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Preface 



Partial differential equations are often used to construct models of the most 
basic theories underlying physics and engineering. For example, the system of 
partial differential equations known as Maxwell's equations can be written on 
the back of a post card, yet from these equations one can derive the entire theory 
of electricity and magnetism, including light. 

Our goal here is to develop the most basic ideas from the theory of partial dif- 
ferential equations, and apply them to the simplest models arising from physics. 
In particular, we will present some of the elegant mathematics that can be used 
to describe the vibrating circular membrane. We will see that the frequencies 
of a circular drum are essentially the eigenvalues from an eigenvector-eigenvalue 
problem for Bessel's equation, an ordinary differential equation which can be 
solved quite nicely using the technique of power series expansions. Thus we 
start our presentation with a review of power series, which the student should 
have seen in a previous calculus course. 

It is not easy to master the theory of partial differential equations. Unlike 
the theory of ordinary differential equations, which relies on the "fundamental 
existence and uniqueness theorem," there is no single theorem which is central 
to the subject. Instead, there are separate theories used for each of the major 
types of partial differential equations that commonly arise. 

However, there are several basic skills which are essential for studying all 
types of partial differential equations. Before reading these notes, students 
should understand how to solve the simplest ordinary differential equations, 
such as the equation of exponential growth dy/dx — ky and the equation of 
simple harmonic motion d 2 y/dx 2 + toy = 0, and how these equations arise in 
modeling population growth and the motion of a weight attached to the ceiling 
by means of a spring. It is remarkable how frequently these basic equations 
arise in applications. Students should also understand how to solve first-order 
linear systems of differential equations with constant coefficients in an arbitrary 
number of unknowns using vectors and matrices with real or complex entries. 
(This topic will be reviewed in the second chapter.) Familiarity is also needed 
with the basics of vector calculus, including the gradient, divergence and curl, 
and the integral theorems which relate them to each other. Finally, one needs 
ability to carry out lengthy calculations with confidence. Needless to say, all 
of these skills are necessary for a thorough understanding of the mathematical 



language that is an essential foundation for the sciences and engineering. 

Moreover, the subject of partial differential equations should not be studied 
in isolation, because much intuition comes from a thorough understanding of 
applications. The individual branches of the subject are concerned with the 
special types of partial differential equations which which are needed to model 
diffusion, wave motion, equilibria of membranes and so forth. The behavior 
of physical systems often suggest theorems which can be proven via rigorous 
mathematics. (This last point, and the book itself, can be best appreciated by 
those who have taken a course in rigorous mathematical proof, such as a course 
in mathematical inquiry, whether at the high school or university level.) 

Moreover, the objects modeled make it clear that there should be a constant 
tension between the discrete and continuous. For example, a vibrating string 
can be regarded profitably as a continuous object, yet if one looks at a fine 
enough scale, the string is made up of molecules, suggesting a discrete model 
with a large number of variables. Moreover, we will see that although a partial 
differential equation provides an elegant continuous model for a vibrating mem- 
brane, the numerical method used to do actual calculations may approximate 
this continuous model with a discrete mechanical system with a large number of 
degrees of freedom. The eigenvalue problem for a differential equation thereby 
becomes approximated by an eigenvalue problem for an n x n matrix where n 
is large, thereby providing a link between the techniques studied in linear al- 
gebra and those of partial differential equations. The reader should be aware 
that there are many cases in which a discrete model may actually provide a 
better description of the phenomenon under study than a continuous one. One 
should also be aware that probabilistic techniques provide an additional compo- 
nent to model building, alongside the partial differential equations and discrete 
mechanical systems with many degrees of freedom described in these pages. 

There is a major dichotomy that runs through the subject — linear versus 
nonlinear. It is actually linear partial differential equations for which the tech- 
nique of linear algebra prove to be so effective. This book is concerned primarly 
with linear partial differential equations yet it is the nonlinear partial differen- 
tial equations that provide the most intriguing questions for research. Nonlinear 
partial differential equations include the Einstein field equations from general 
relativity and the Navier-Stokes equations which describe fluid motion. We hope 
the linear theory presented here will whet the student's appetite for studying 
the deeper waters of the nonlinear theory. 

The author would appreciate comments that may help improve the next 
version of this short book. He hopes to make a list of corrections available at 
the web site: 

http : //www . math . ucsb . edu/~moore 

Doug Moore 
March, 2003 
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Chapter 1 

Power Series 



1.1 What is a power series? 

Functions are often represented efficiently by means of infinite series. Examples 
we have seen in calculus include the exponential function 

11 °° 1 

e* = l + x+-x>+-x 3 + ... = J2-iX n , (1-1) 

as well as the trigonometric functions 



cosx = 1 - i* 2 + ix 4 - • • • = JL,- 



fe=0 



2fe 

(2fc)f 



and 

1 , 1 



1 1 

= * - r 3 + y* 5 — ■ = E(-D fc j^Tiy- 



2k+l 



An infinite series of this type is called a power series. To be precise, a power 
series centered at x n is an infinite sum of the form 

oo 

a + ai(x - x ) + a 2 (x - x ) 2 H = ^ a n( x ~ x o) n , 

where the o„'s are constants. In advanced treatments of calculus, these power 
series representations are often used to define the exponential and trigonometric 
functions. 

Power series can also be used to construct tables of values for these functions. 
For example, using a calculator or PC with suitable software installed (such as 
Mathematica), we could calculate 

1 + 1 + ^ l2 = E^ in = 2 - 5 > 

n=0 



1 



4 

1 + 1 + ^l 2 + il 3 + = Y -X = 2.70833, 
2! 3! 4! ^ n! 

8 12 

V — 1" = 2.71806, V — 1" = 2.71828, 
^ n! ^ n! 

n=0 n=0 

As the number of terms increases, the sum approaches the familiar value of the 
exponential function e x at x = 1. 

For a power series to be useful, the infinite sum must actually add up to a 
finite number, as in this example, for at least some values of the variable x. We 
let sjv denote the sum of the first (N + 1) terms in the power series, 

N 

sn = ao + ai(x - x ) + a 2 (x - x ) 2 H h a N (x - x ) N = ^ a n (x - x ) n , 

and say that the power series 



y a n (x - x ) r 



n=0 



converges if the finite sum s n gets closer and closer to some (finite) number as 
N — > oo. 

Let us consider, for example, one of the most important power series of 
applied mathematics, the geometric series 



oo 

x n . 

n=0 



i + x + x 2 + x 3 h — = y~] 

In this case we have 

s N = 1 + x + x 2 + x 3 H h/, xs N = x + x 2 + x 3 + x 4 h x N+1 , 

N+l 1 x 



S^v — XSjy = I — X , S]v 



1-a; 

If | a; | < 1, then x N+1 gets smaller and smaller as N approaches infinity, and 
hence 

lim x N+1 = 0. 
Substituting into the expression for s?f, we find that 

lim sn 



N^oo 1 — X 

Thus if |x| < 1, we say that the geometric series converges, and write 



oo 1 

5>" = T^ 



n=0 



2 



On the other hand, if \x\ > 1, then x N+1 gets larger and larger as N ap- 
proaches infinity, so liniAr^oo x N+1 does not exist as a finite number, and neither 
does liniAT^oo s^. In this case, we say that the geometric series diverges. In 
summary, the geometric series 

x i 

^ x n converges to when \x\ < 1, 

n=o x 

and diverges when |x| > 1. 

This behaviour, convergence for \x\ < some number, and divergences for 
\x\ > that number, is typical of power series: 

Theorem. For any power series 

oo 

a + ai(x - x ) + a 2 (x - x ) 2 H = ^ a n( x - x aT, 

n=0 

there exists R, which is a nonnegative real number or oo, such that 

1. the power series converges when \x — x \ < R, 

2. and the power series diverges when \x — x n \ > R. 



We call R the radius of convergence . A proof of this theorem is given in more 
advanced courses on real or complex analysis. 1 
We have seen that the geometric series 



1 + x + x 2 + x 3 + 



n=0 



has radius of convergence R = 1. More generally, if b is a positive constant, the 
power series 



1 + 



+ ■ 



£(!)" 



(1.2) 



has radius of convergence b. To see this, we make the substitution y — x/b, 
and the power series becomes J^^Lo^™' w hich we already know converges for 
\y\ < 1 and diverges for \y\ > 1. But 



\y\ < 1 

\y\ > 1 ^ 



< 1 ^ \x\ < b, 



> 1 



Id > b. 



1 Good references for the theory behind convergence of power series are Edward D. 
Gaughan, Introduction to analysis, Brooks/Cole Publishing Company, Pacific Grove, 1998 
and Walter Rudin, Principles of mathematical analysis, third edition, McGraw-Hill, New 
York, 1976. 
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Thus for | a; | < b the power series (1.2) converges to 

1 1 b 



1 — y 1 — (x/b) b — x' 

while for |x| > b, it diverges. 

There is a simple criterion that often enables one to determine the radius of 
convergence of a power series. 

Ratio Test. The radius of convergence of the power series 

oo 

a + a 1 (x - x ) + a 2 (x - x ) 2 H = ^ a n (x - x ) n 

is given by the formula 

R= lim '"" 



n^oo \a n+1 \ 

so long as this limit exists. 

Let us check that the ratio test gives the right answer for the radius of conver- 
gence of the power series (1.2). In this case, we have 

1_ K| 1/6" _ b n+1 _ u 

a "-6«' S ° |a n+i r V&" +1 " 

and the formula from the ratio test tells us that the radius of convergence is 
R = b, in agreement with our earlier determination. 
In the case of the power series for e x , 



L — ' n 



n=0 

in which a n = 1/n!, we have 

Kj 1/n! = (n+l)\ 

\a n+1 \ l/(n + l)! n! 

and hence 



= n + 1, 



i? = lim a " = lim (n + 1) = oo, 

n^oo |a n _|_i| n^oo 

so the radius of convergence is infinity. In this case the power series converges 
for all x. In fact, we could use the power series expansion for e x to calculate e x 
for any choice of x. 

On the other hand, in the case of the power series 



n=0 



4 



in which a„ — n\, we have 



R = lim -r—^-r = lim = 0. 



|cin+i| (n + 1)! n+1' n ^°° \a n +i\ n^oo^n+1 

In this case, the radius of convergence is zero, and the power series does not 
converge for any nonzero x. 

The ratio test doesn't always work because the limit may not exist, but 
sometimes one can use it in conjunction with the 

Comparison Test. Suppose that the power series 

oo oo 
n=0 n=0 

have radius of convergence R\ and R2 respectively. If \a n \ < \b n \ for all n, then 
Ri > Ri- If K| > \K\ for all n, then R 1 < R 2 . 

In short, power series with smaller coefficients have larger radius of convergence. 
Consider for example the power series expansion for cos x, 



1 

1 + to - -x 2 + Ox 3 + -x 4 - • • • = — , 



k=0 



2k 

{2k)V 



In this case the coefficient a n is zero when n is odd, while a n = ±l/n! when 
n is even. In either case, we have \a n \ < l/n\. Thus we can compare with the 
power series 

n=0 

which represents e x and has infinite radius of convergence. It follows from the 
comparison test that the radius of convergence of 



00 



x 2k 



must be at least large as that of the power series for e x , and hence must also be 
infinite. 

Power series with positive radius of convergence are so important that there 
is a special term for describing functions which can be represented by such power 
series. A function f(x) is said to be real analytic at xq if there is a power series 



a n(x - x ) n 



n=Q 

about x with positive radius of convergence R such that 



f(x) = a n (x - x ) n , for \x - x \ < R. 



n=0 
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For example, the functions e x is real analytic at any Xq. To see this, we 
utilize the law of exponents to write e x — e x o e x-x and apply (1.1) with x 
replaced by x — x : 

~ 1 ~ e *o 

e x = e x ° } — (x - x ) n = )a n (x - x ) n , where a n = — . 

-' TV. ^— ' TV. 

This is a power series expansion of e x about x with infinite radius of con- 
vergence. Similarly, the monomial function f(x) = x n is real analytic at x 
because 

x n = (x - x + x Q ) n = V U ' ... xl~\x - x a y 
^ inn — 1)\ 

i=0 v ' 

by the binomial theorem, a power series about x a in which all but finitely many 
of the coefficients are zero. 

In more advanced courses, one studies criteria under which functions are 
real analytic. For our purposes, it is sufficient to be aware of the following facts: 
The sum and product of real analytic functions is real analytic. It follows from 
this that any polynomial 

P(x) = a + a\X + a 2 x 2 + • • • + a n x n 

is analytic at any x . The quotient of two polynomials with no common factors, 
P(x)/Q(x), is analytic at xq if and only if x n is not a zero of the denominator 
Q(x). Thus for example, l/(x — 1) is analytic whenever x ^ 1, but fails to be 
analytic at xq = 1. 

Exercises: 

1.1.1. Use the ratio test to find the radius of convergence of the following power 
series: 



a. b. 

OO - oo _. 

c. Y (z-2)", d. V— {x-w) n , 

oo oo 

e. £(7* -14)", f. E^ 3 *- 6 )"- 

n=0 n=0 

1.1.2. Use the comparison test to find an estimate for the radius of convergence 
of each of the following power series: 



.. E^-, b. 



2fc 



fc=0 fe=0 



G 



1.1.3. Use the comparison test and the ratio test to find the radius of convergence 
of the power series 

^ [ ' (ml) 2 \2J 

1.1.4. Determine the values of xq at which the following functions fail to be real 
analytic: 

1 x 
a. fW = JZ7Z> b - 9(x) = ^j, 

4 1 

c. h(x) — —. -, d. 6(x) — —5 — = — 

K ' x 4 - 3x 2 + 2 ' Yy ' x 3 - hx 2 + 6x 

1.2 Solving differential equations by means of 
power series 

Our main goal in this chapter is to study how to determine solutions to differ- 
ential equations by means of power series. As an example, we consider our old 
friend, the equation of simple harmonic motion 

which we have already learned how to solve by other methods. Suppose for the 
moment that we don't know the general solution and want to find it by means 
of power series. We could start by assuming that 

oo 

y = a + aix + a 2 x 2 + a 3 x 3 H = a n x n . (1.4) 

It can be shown that the standard technique for differentiating polynomials term 
by term also works for power series, so we expect that 

— = a\ + 2a 2 x + 3a 3 x 2 + • • • = > na n x n ~ l . 
ax ^— ; 

n—l 

(Note that the last summation only goes from 1 to oo , since the term with n = 
drops out of the sum.) Differentiating again yields 

d 2 °° 

^§ = 2a 2 + 3 • 2a 3 x + 4 • 3a 4 x 2 H = ^ n(n - l)a„x n " 2 . 

71=2 



We can replace n by m + 2 in the last summation so that 

,o OO OO 

dx 2 



72 

"' V = (m + 2)(m + 2-l)a m+2 x m+2 - 2 = Y,(m + 2)(m + l)a m+2 x r - 



iyi+2=2 m=0 
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The index m is a "dummy variable" in the summation and can be replaced by 
any other letter. Thus we are free to replace m by n and obtain the formula 

0=^(n + 2)(n + lK +2 x". 

n=0 

Substitution into equation(1.3) yields 

00 00 
^2(n + 2)(n + l)a n+2 x n + 

n=0 n=0 

or 

00 

Y^[(n + 2)(n + l)a„+ 2 + ajz™ = 0. 

n=0 

Recall that a polynomial is zero only if all its coefficients are zero. Similarly, a 
power series can be zero only if all of its coefficients are zero. It follows that 

(n + 2)(n + l)a n+2 + a n = 0, 



or 



l n+2 



(n + 2)(n + l)' 



(1.5) 



This is called a recursion formula for the coefficients o n . 

The first two coefficients ao and a\ in the power series can be determined 
from the initial conditions, 

2/(0) = a , £(0) = ai . 

Then the recursion formula can be used to determine the remaining coefficients 
by the process of induction. Indeed it follows from (1.5) with n — that 

a 1 

08 = -— = --O0. 

Similarly, it follows from (1.5) with n = 1 that 



and with n = 2 that 



ai I 

a 3 = -^ = -^«l, 



a 2 11 1 



Continuing in this manner, we find that 



_ (-1)" _ ("!)" 

a2n ~ W a °' a2n+1 (2^TT)! ai - 
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Substitution into (1.4) yields 

y = a + aix - ^a x 2 - Ty^x 3 + ^j ^ 4 H 

We recognize that the expressions within parentheses are power series expan- 
sions of the functions sin x and cos x, and hence we obtain the familiar expression 
for the solution to the equation of simple harmonic motion, 

y = a cos x + a\ sin x. 

The method we have described — assuming a solution to the differential equa- 
tion of the form 

oo 

y(x) = a ™ x ™ 

n=0 

and solve for the coefficients a n — is surprisingly effective, particularly for the 
class of equations called second-order linear differential equations. 

It is proven in books on differential equations that if P(x) and Q(x) are well- 
behaved functions, then the solutions to the "homogeneous linear differential 
equation" 

can be organized into a two-parameter family 

y = a y (x) +a 1 y 1 (x), 

called the general solution. Here yo(x) and y\{x) are any two nonzero solutions, 
neither of which is a constant multiple of the other. In the terminology used in 
linear algebra, we say that they are linearly independent solutions. As ao and 
ai range over all constants, y ranges throughout a "linear space" of solutions. 
We say that yo{x) and yi(x) form a basis for the space of solutions. 

In the special case where the functions P(x) and Q(x) are real analytic, 
the solutions yo(x) and y\(x) will also be real analytic. This is the content of 
the following theorem, which is proven in more advanced books on differential 
equations: 

Theorem. If the functions P(x) and Q(x) can be represented by power series 

oo oo 

P(x) = ^2p n (x- x ) n , Q(x) = ^2q n (x- x ) n 

n=0 n=0 

with positive radii of convergence R\ and R2 respectively, then any solution 
y(x) to the linear differential equation 

g + P W f + Q(x)„ = 



9 



can be represented by a power series 

oo 
n=0 

whose radius of convergence is > the smallest of R x and R 2 . 

This theorem is used to justify the solution of many well-known differential 
equations by means of the power series method. 

Example. Hermite's differential equation is 

where p is a parameter. It turns out that this equation is very useful for treating 
the simple harmonic oscillator in quantum mechanics, but for the moment, we 
can regard it as merely an example of an equation to which the previous theorem 
applies. Indeed, 

P(x) = -2x, Q(x) = 2p, 

both functions being polynomials, hence power series about Xq = with infinite 
radius of convergence. 

As in the case of the equation of simple harmonic motion, we write 



y = ^a n x n . 



n=0 

We differentiate term by term as before, and obtain 
d V \^ .„-i d 2 y 



= ^2na n x n 1 , — = ^2n(n- l)a n x" 



,n-2 

ax < dx z 

n=l n=2 



Once again, we can replace n by to + 2 in the last summation so that 
dx 2 



oo 



j2 00 

^2 (m + 2)(to + 2 - l)a m+2 x m+2 - 2 = ^ (m + 2)(m + l)a m+2 x n 



m+2=2 m=0 

and then replace m by n once again, so that 
d 2 y 



Note that 



= £(n+2)(n + l)a„ + 2S n . (1.7) 

n=0 



A 

-2x4- = V -2na n x n , (1.8) 
dx ^ 

n=0 
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while 

CO 

2py = ^2pa n x n . (1.9) 
Adding together (1.7), (1.8) and (1.9), we obtain 

_ 2x -^ + 2py = J2( n + 2)(n + l)a n+2 x n + ^(-2n + 2p)a n x n . 
If y satisfies Hermite's equation, we must have 

oo 

= ^2[{n + 2)(n + l)a„ +2 (-2n + 2p)a„]x™. 

n=0 

Since the right-hand side is zero for all choices of x, each coefficient must be 
zero, so 

(n + 2)(n + l)a n+2 + (-2n + 2p)a n = 0, 
and we obtain the recursion formula for the coefficients of the power series: 

2n-2p 

°"+ 2 = (n+2)(n+l) a - (L10) 

Just as in the case of the equation of simple harmonic motion, the first two 
coefficients cto and a\ in the power series can be determined from the initial 
conditions, 

1/(0) = oo, £(0) = ai. 

The recursion formula can be used to determine the remaining coefficients in 
the power series. Indeed it follows from (1.10) with n = that 

2p 

Similarly, it follows from (1.10) with n = 1 that 



2 - 2p 2(p - 1) 

"3 = 

and with n = 2 that 



a3 = ^2~ ai = 3!" 



4 - 2p 2(2 - p) -2p 2 2 p(p - 2) 

Continuing in this manner, we find that 

6~2p 2(3 -p) 2(1 -p) 2 2 (p-l)(p-3) 
° 5 = T^ 3 = —A 3i~ 01 = 5! ° 1 ' 
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8-2p 2(3-p)2 2 (p-2)p 

a 4 — — 7, — r 7i a 



2 3 p(p-2)(p-4) 



6-5-2 * 6-5 
and so forth. Thus we find that 



4! 



6! 



a , 



y = a a 



2p 2 , 2 2 p(p- 2) 4 2 3 p(p-2)( P -4) 6 
2! + 4! 6! + 



+(/| g _ g^l) 8 , + 2'( P - 1)^-3)^ 

2 3 (p-l)(p-3)(p-5) 7 

7! + " 
We can now write the general solution to Hcrmitc's equation in the form 

y = a y {x) +a 1 y 1 (x), 

where 



Vo{x) 



2p 2 2 p(p-2) _ 2 3 p(p-2)(p-4) 

2! 4! 6! + "' 



and 



3! 



5! 



7! 



For a given choice of the parameter p, we could use the power series to construct 
tables of values for the functions yo(x) and yi(x). Tables of values for these 
functions are found in many "handbooks of mathematical functions." In the 
language of linear algebra, we say that yo(x) and yi(x) form a basis for the 
space of solutions to Hermite's equation. 

When p is a positive integer, one of the two power series will collapse, yielding 
a polynomial solution to Hcrmitc's equation. These polynomial solutions are 
known as Hermite polynomials. 



Another Example. Legendre's differential equation is 



dx 2 



2x^+p(p+l)y = 0, 



(1.11) 



where p is a parameter. This equation is very useful for treating spherically 
symmetric potentials in the theories of Newtonian gravitation and in electricity 
and magnetism. 

To apply our theorem, we need to divide by 1 — x 2 to obtain 



Thus we have 



dx 2 



P(x) 



2x dy p(p+l) 

, _ ,.2 



1 — x 2 dx 
2x 



l-x 



2 : 



Q(x) 



p(p+i) 

l-x 2 ' 
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Now from the preceding section, we know that the power series 

1 + u + u 2 + u 3 + ■ ■ ■ converges to — - — 

1 — u 

for |tt| < 1. If we substitute u = x 2 , we can conclude that 

- = 1 + x 2 + x 4 + x 6 + ■ ■ ■ , 



1 - 

the power series converging when |x| < 1. It follows quickly that 

2x 

P(x) = = -2x - 2x 3 -2x 5 

1 — x z 



and 

l-x 



Q( X ) = T _ 2 =P(P+1)+P(P+ 1)X 2 + p(p +1)X 4 + 



Both of these functions have power series expansions about xq = which con- 
verge for \x\ < 1. Hence our theorem implies that any solution to Lcgcndre's 
equation will be expressible as a power series about x$ = which converges for 
\x\ < 1. However, we might suspect that the solutions to Legendre's equation 
to exhibit some unpleasant behaviour near x = ±1. Experimentation with nu- 
merical solutions to Legendre's equation would show that these suspicions are 
justified — solutions to Legendre's equation will usually blow up as x — > ±1. 

Indeed, it can be shown that when p is an integer, Legendre's differential 
equation has a nonzero polynomial solution which is well-behaved for all x, but 
solutions which are not constant multiples of these Legendre polynomials blow 
up as x — * ±1. 

Exercises: 

1.2.1. We would like to use the power series method to find the general solution 
to the differential equation 

^i-44 + 12, = 0, 
ax z ax 

which is very similar to Hcrmite's equation. So we assume the solution is of the 
form 

oo 

V = ^2 a n x n , 

n=0 

a power series centered at 0, and determine the coefficients a n . 

a. As a first step, find the recursion formula for a n+2 in terms of a n . 

b. The coefficients ao and a\ will be determined by the initial conditions. Use 
the recursion formula to determine a n in terms of a and a\, for 2 < n < 9. 

c. Find a nonzero polynomial solution to this differential equation. 
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d. Find a basis for the space of solutions to the equation. 

e. Find the solution to the initial value problem 

g-*£ + i*-o, ,(0,-0, *(0>-l. 

f. To solve the differential equation 

d2 y M ?\ d y , 10 n 
— - 4(x - 3)— + 12y = 0, 

dx z ax 

it would be most natural to assume that the solution has the form 

00 

y = Y / a n (x-3) n . 

n=0 

Use this idea to find a polynomial solution to the differential equation 

1.2.2. We want to use the power series method to find the general solution to 
Legendre's differential equation 

Once again our approach is to assume our solution is a power series centered at 
and determine the coefficients in this power series. 

a. As a first step, find the recursion formula for a n+ 2 in terms of a n . 

b. Use the recursion formula to determine a n in terms of ciq and ai, for 2 < n < 
9. 

c. Find a nonzero polynomial solution to this differential equation, in the case 
where p = 3. 

d. Find a basis for the space of solutions to the differential equation 

)^L-2x d ^- 
dx 2 dx 



1.2.3. The differential equation 

2 d 2 y dy 
{l - x) dx^ X Tx +PV = ^ 



where p is a constant, is known as Chebyshev's equation. It can be rewritten in 
the form 

d 2 y dy x p 2 



P{x)-f+Q{x)y = 0, where P{x) = -- Q(x) 



dx 2 dx ' 1 — x 2 1 — x 2 
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a. If P(x) and Q(x) are represented as power series about xq = 0, what is the 
radius of convergence of these power series? 

b. Assuming a power series centered at 0, find the recursion formula for a rl +2 
in terms of a n . 

c. Use the recursion formula to determine a n in terms of a and ai, for 2 < n < 
9. 

d. In the special case where p — 3, find a nonzero polynomial solution to this 
differential equation. 

e. Find a basis for the space of solutions to 
1.2.4. The differential equation 

■^ + ^jz = Xz (1.12) 

arises when treating the quantum mechanics of simple harmonic motion. 

a. Show that making the substitution z = er x / 2 y transforms this equation into 
Hermite's differential equation 

£-*£ + (*-i).-o. 

b. Show that if A = 2n+l where n is a nonnegative integer, (1.12) has a solution 
of the form z = e~ x / 2 P n {x), where P n (x) is a polynomial. 

1.3 Singular points 

Our ultimate goal is to give a mathematical description of the vibrations of a 
circular drum. For this, we will need to solve Bessel's equation, a second-order 
homogeneous linear differential equation with a "singular point" at 0. 
A point x is called an ordinary point for the differential equation 

g + P (a5 )| + 0(^ = (1.13) 

if the coefficients P(x) or Q(x) are both real analytic at x — x , or equivalently, 
both P(x) or Q(x) have power series expansions about x — x n with positive 
radius of convergence. In the opposite case, we say that x is a singular point; 
thus xo is a singular point if at least one of the coefficients P(x) or Q(x) fails 
to be real analytic at x — x - A singular point is said to be regular if 

(x — x a )P(x) and (x — x ) 2 Q(x) 
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are real analytic. 

For example, xq = 1 is a singular point for Lcgendrc's equation 

d 2 y 2x dy | p(p+l) _ ^ 
dx 2 1 — x 2 dx I — x 2 

because 1 — x 2 — > as x — > 1 and hence the quotients 

2, md p(p + l) 



1 — x 2 1 — x 2 

blow up as a; — > 1 , but it is a regular singular point because 

-2x 2x 



(x-l)P(x) = (x-l) 



1 — x 2 x + 1 



and 



(x-l) 2 Q(x) = (x-l) 2PiP+1) _P(P+1)(1-^) 



1 — X 2 1 + X 

are both real analytic at Xo = 1- 

The point of these definitions is that in the case where x = xo is a regular 
singular point, a modification of the power series method can still be used to 
find solutions. 

Theorem of Frobenius. If xo is a regular singular point for the differential 
equation 

d 2 y | P f x \ d v 

dx 2 dx 

then this differential equation has at least one nonzero solution of the form 



,.2 +P(x)± + Q(x)y = 0, 



y(x) = (x-x ) r ^2a n (x-x a ) n , (1.14) 



n=0 

where r is a constant, which may be complex. If (x — x )P(x) and (x~xo) 2 Q(x) 
have power series which converge for \x — x \ < R then the power series 

oo 

^2 a n (x - x ) n 

n=0 

will also converge for \x — xq\ < R. 

We will call a solution of the form (1.14) a generalized power series solution. 
Unfortunately, the theorem guarantees only one generalized power series solu- 
tion, not a basis. In fortuitous cases, one can find a basis of generalized power 
series solutions, but not always. The method of finding generalized power series 
solutions to (1.13) in the case of regular singular points is called the Frobenius 
method. 2 



2 For more discussion of the Frobenius method as well as many of the other techniques 
touched upon in this chapter we refer the reader to George F. Simmons, Differential equations 
■with applications and historical notes, second edition, McGraw-Hill, New York, 1991. 
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The simplest differential equation to which the Theorem of Frobcnius applies 
is the Cauchy-Euler equidimensional equation. This is the special case of (1.13) 
for which 

P{x) = P -, Q{x) = \, 
x x z 

where p and q are constants. Note that 

xP(x) = p and x 2 Q(x) = q 

are real analytic, so x — is a regular singular point for the Cauchy-Euler 
equation as long as either p or q is nonzero. 

The Frobcnius method is quite simple in the case of Cauchy-Euler equations. 
Indeed, in this case, we can simply take y(x) = x r , substitute into the equation 
and solve for r. Often there will be two linearly independent solutions y\{x) = 
x Tl and 2/2(2:) = x r2 of this special form. In this case, the general solution is 
given by the superposition principle as 

y = cix Tl + c 2 x T ' 2 . 

For example, to solve the differential equation 

2 d 2 y , dy 
aar ax 

we set y — x r and differentiate to show that 

dy/dx = rx r ~ 1 => x{dy/dx) — rx r , 

d 2 y/dx 2 = r(r — l)x r ~ 2 => x 2 (d 2 y/dx 2 ) = r(r - l)x r . 

Substitution into the differential equation yields 

r(r-l)x r +Arx r + 2x r = 0, 

and dividing by x r yields 

r(r - 1) + 4r + 2 = or r 2 + 3r + 2 = 0. 

The roots to this equation are r = — 1 and r = —2, so the general solution to 
the differential equation is 

-1 , -2 C l , C 2 

y = c\x + c 2 x = h -^r. 

x x l 

Note that the solutions y\{x) = x^ 1 and y 2 {x) = x~ 2 can be rewritten in the 
form 

00 00 
2/i (x) = x^ 1 ^2 a nX n , y 2 (x) = x~ 2 ^2b n x n , 

n=0 n=0 

where ao = 60 = 1 and all the other a„'s and 6„'s are zero, so both of these 
solutions are generalized power series solutions. 
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On the other hand, if this method is applied to the differential equation 

2 d2 V . o d V , n 

x d^ + ix T x + y = ^ 

we obtain 

r(r - 1) + 3r + 1 = r 2 + 2r + 1, 

which has a repeated root. In this case, we obtain only a one-parameter family 
of solutions 

y = cx^ 1 . 

Fortunately, there is a trick that enables us to handle this situation, the so-called 
method of variation of parameters. In this context, we replace the parameter c 
by a variable v(x) and write 

y = v(x)x^ 1 . 

Then 

^ = v'(x)x~ 1 - v(x)x~ 2 , = v"(x)x~ 1 - 2v'(x)x~ 2 + 2v(x)x~ 3 . 

dx dx z 

Substitution into the differential equation yields 

x 2 (v"(x)x~ 1 - 2v'(x)x- 2 + 2v(x)x~ 3 ) + 3x(v'(x)x- 1 - v(x)x~ 2 ) + v(x)x~ 1 = 0, 
which quickly simplifies to yield 

xv"(x) + v'(x) = 0, — - = --, loglt/l = -log Id +a, v'=—, 
v' x x 

where a and c 2 are constants of integration. A further integration yields 

v = c 2 log |a;j +ci, so y = (c 2 log \x\ + ci)a; _1 , 
and we obtain the general solution 

1 log | a; | 
y = ci- +c 2 - 



x x 

In this case, only one of the basis elements in the general solution is a generalized 
power series. 

For equations which are not of Cauchy-Euler form the Frobenius method is 
more involved. Let us consider the example 

which can be rewritten as 

^y + P{x) ^ + Q{x)y = , whcrc p(x) = ± QW = Yx- 
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One easily checks that x = is a regular singular point. We begin the Frobenius 
method by assuming that the solution has the form 



CO 

n+r 



y = x r a « x " = a " 

Then 



.x 

n=0 n=0 



i oo ,2 

-£ = ^(n + r)a n x n+r -\ || = ]T(n + r)(n + r - l)a n x n 



„n+r-2 

n=0 ™~ n=0 

and 



j2 00 

2x^ = ^ 2(n + r)(n + r - l)a n x n+r - 1 . 
Substitution into the differential equation yields 

oo oo oo 

2(n + r)(n + r - l)o n a; n+r - 1 + ^(n + v)a n x n+r - x + £ a„a;" +r = 0, 

n— n— n— 

which simplifies to 



J^(2n + 2r - l)(n + r)a„a;"" 1 + ^ a n x n 

.n=0 n=0 



= 0. 



We can divide by x r , and separate out the first term from the first summation, 
obtaining 

(2r - l)ra z _1 + J^(2n + 2r - l)(n + r)a„or n_1 + a„a;" = 0. 

n— 1 n— 

If we let n = m + 1 in the first infinite sum, this becomes 

oo oo 

(2r - ^raox" 1 + ^ (2m + 2r + l)(m + r + l)a m+1 x m + ^ a « x ™ = °- 

m — n— 

Finally, we replace m by n, obtaining 

oo oo 

(2r - ^raoa;- 1 + ^(2n + 2r + l)(n + r + l)a n+ iic n + ^ a„a;" = 0. 

n=0 n=0 

The coefficient of each power of x must be zero. In particular, we must have 

(2r-l)rao = 0, (2n + 2r + l)(n + r + l)a n+1 + a„ = 0. (1.16) 

If ao = 0, then all the coefficients must be zero from the second of these equa- 
tions, and we don't get a nonzero solution. So we must have a ^ and hence 

(2r - l)r = 0. 
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This is called the indicial equation. In this case, it has two roots 



ri = 0, r 2 = -. 

The second half of (1.16) yields the recursion formula 

1 



(2n + 2r + l)(n + r + 1) 



for n > 0. 



We can try to find a generalized power series solution for either root of the 
indicial equation. If r = 0, the recursion formula becomes 

1 

fl n+l — - 



(2n + l)(n + l) 



Given a = 1, we find that 

<H = -1, 



02 



3-2 ai " 3-2' 



1 



1 



«3 



5-3 02 (5-3)(3-2)' 
and so forth. In general, we would have 

a n = (-1)"- 



04 



7-4° 3 (7- 5- 3)4! : 



(2n-l)(2n-3)---l-n!' 
One of the generalized power series solution to (1.15) is 



y\{x) = x° 



l-x + 
l-x + 



1 



1 



3-2 (5-3)(3!) (7-5-3)4! 

1 o 1 , . 1 



3-2 (5-3)(3!) (7-5-3)4! 
If r = 1/2, the recursion formula becomes 

1 1 



(2n+2)(n+(l/2) + l) 
Given a = 1, we find that 



(n+l)(2n + 3) 



On. 



«i = 



fl 2 = — ^— -ai — 



1 



2-5 i_ 2-5-3' 
1 



«3 = -7T^ a 2 



and in general, 



On = (-1)" 



3-7 ' 3! -(7-5-3)' 
1 

n!(2n+l)(2n-l)---l-n!' 
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We thus obtain a second generalized power series solution to (1.15): 



y 2 (x) = x 1 ' 2 



1 



1 



1 



3 2-5-3 3! • (7 • 5 • 3) 
The general solution to (1.15) is a superposition of y\(x) and yzix): 

1 o 1 1 



y = c-i 



1 — x + -——x" 
3 • 2 



(5-3)(3!)' 



(7-5-3)4 



1 



1 



2-5-3 3! • (7 • 5 • 3) 

We obtained two linearly independent generalized power series solutions in 
this case, but this does not always happen. If the roots of the indicial equation 
differ by an integer, we may obtain only one generalized power series solution. 
In that case, a second independent solution can then be found by variation of 
parameters, just as we saw in the case of the Cauchy-Euler cquidimcnsional 
equation. 

Exercises: 

1.3.1. For each of the following differential equations, determine whether x = 
is ordinary or singular. If it is singular, determine whether it is regular or not. 

a. y" + xy' + (1 - x 2 )y = 0. 

b. y" + (l/x)y' + (1 - (l/x 2 ))y = 0. 

c. x 2 y" + 2xy' + (cos x)y = 0. 

d. x 3 y" + 2xy' + (cosx)y = 0. 

1.3.2. Find the general solution to each of the following Cauchy-Euler equations: 

a. x 2 d 2 y/dx 2 - 2xdy/dx + 2y = 0. 

b. x 2 d 2 y/dx 2 - xdy/dx + y = 0. 

c. x 2 d 2 y/dx 2 - xdy/dx + lOy = 0. 
(Hint: Use the formula 



a+bi _ a 



a(„logx\bi _ a ib log x _ a 



= x a x bl = x a {e losx r = x a e 



= x a [cos(6 log x) + i sin(6 log x) 



to simplify the answer.) 



1.3.3. We want to find generalized power series solutions to the differential 
equation 

Q d 2 y dy 

3x dx- 2 + d-x +y = ° 
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by the method of Frobcnius. Our procedure is to find solutions of the form 

oo oo 

y = x r a n x n = Y a n x n+r , 
where r and the a„'s are constants. 

a. Determine the indicial equation and the recursion formula. 

b. Find two linearly independent generalized power series solutions. 

1.3.4. To find generalized power series solutions to the differential equation 

„ d 2 y dy 

2x t4 + -f + X V = 

dx z ax 

by the method of Frobenius, we assume the solution has the form 

oo 

y = J2a n x n+r , 

n=0 

where r and the a„'s are constants. 

a. Determine the indicial equation and the recursion formula. 

b. Find two linearly independent generalized power series solutions. 

1.4 Bessel's differential equation 

Our next goal is to apply the Frobenius method to Bessel's equation, 

^64W-p 2 )</ = 0, (1.17) 



dx \ dx 

an equation which is needed to analyze the vibrations of a circular drum, as we 
mentioned before. Here p is a parameter, which will be a nonnegative integer 
in the vibrating drum problem. Using the Leibniz rule for differentiating a 
product, we can rewrite Bessel's equation in the form 

or equivalently as 

g + P(x)| + «.)-«. 

where 

P(x) = - and Q[x)=^-^L. 
x x z 

Since 

xP(x) = 1 and x 2 Q(x) = x 2 — p 2 , 
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wc sec that x = is a regular singular point, so the Frobenius theorem implies 
that there exists a nonzero generalized power series solution to (1.17). 

To find such a solution, we start as in the previous section by assuming that 



y = J2a n x n+r . 



n=0 

Then 

2 = + rKa^- 1 , = + r)a n z"+' 

n=0 n=0 



dy 

dx \~ dx 



r— 1 =J2(n + r) 2 a n x n+r -\ 



and thus 



x— a; 



= ^(n + r) 2 a„a;"+ r . (1.18) 



71=0 



On the other hand, 

00 00 

X 2 y = a nX n+r+2 = E a m-2X m+T , 

n— m— 2 

where we have set m = n + 2. Replacing m by n then yields 



x 2 



y = J2^-2X n+r . (1.19) 



n=2 

Finally, we have, 



-P 2 y=-Y.P 2a ^ n+r - ( L2 °) 

Adding up (1.18), (1.19), and (1.20), we find that if y is a solution to (1.17), 

oo oo oo 

£(n + r) 2 fl „ I "+ r + ^ a„_ 2 x"+'' - £ p 2 a n x n+r = 0. 

n— n— 2 n— 

This simplifies to yield 

oo oo 

£ [(n + rf - p 2 ]a n x n+r + £ a„_ 2 .*"+'' = 0, 
or after division by x r , 

oo oo 

^[(n + r) 2 - p 2 ]a„a;" + ^ a„_ 2 a;" = 0. 

n=0 n=2 
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Thus we find that 

oo 

(r 2 - p 2 )a + [{r + l) 2 - p 2 }a lX + ^{[(n + rf - p 2 }a n + a n ^}x n = 0. 

n=2 

The coefficient of each power of x must be zero, so 

(r 2 -p 2 )a = 0, [(r + l) 2 -p 2 ] ai = 0, [{n+rf -p 2 ]a n + a n - 2 = for n > 2. 
Since we want oq to be nonzero, r must satisfy the indicial equation 

(r 2 -p 2 ) = 0, 

which implies that r — ±p. Let us assume without loss of generality that p > 
and take r = p. Then 

[0 + l) 2 - p 2 \a x = =4> (2p + l)oi = => ai = 0. 

Finally, 

[(n + p) 2 -p 2 ]a„ + a„_ 2 = [n 2 + 2np]a„ + a„_ 2 = 0, 

which yields the recursion formula 

a« = -o ~, — o a "-2- ( L21 ) 

2np + n" 1 

The recursion formula implies that a n = if n is odd. 

In the special case where p is a nonnegative integer, we will get a genuine 
power series solution to Bessel's equation (1.17). Let us focus now on this 
important case. If we set 

1 

a °- 2^!' 



we obtain 

-go _J L r-i^lV +2 

° 2 4p + 4 4(p+l)2Pp! 1 >\2J l!(p+l)!' 

« / 1 \ P+2 * 

-a 2 1 /1\ 1 

CI4 



16 8(p+2) \2J l!(p+ 1)! 

(-1) 2 ' ^ ' 



2(p+2) V2 / l!(p+l)! v ' \2j 2!(p + 2)! : 
and so forth. The general term is 

p+2m 



(-, x p+zm .J 
2 J m!(p + m)! 
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Figure 1.1: Graph of the Bessel function Jq(x). 



Thus we finally obtain the power series solution 



y = 



(!)'£<-» 



m 



m\{p+m)\ V2 




m=0 



The function defined by the power series on the right-hand side is called the 
p-th Bessel function of the first kind, and is denoted by the symbol J p (x). For 
example, 



Using the comparison and ratio tests, we can show that the power series ex- 
pansion for J p {x) has infinite radius of convergence. Thus when p is an integer, 
Bessel's equation has a nonzero solution which is real analytic at x = 0. 

Bessel functions are so important that Mathematica includes them in its 
library of built-in functions. 3 Mathematica represents the Bessel functions of 
the first kind symbolically by BesselJ[n,x] . Thus to plot the Bessel function 
J n (x) on the interval [0, 15] one simply types in 

n=0; Plot [ Bessel J [n,x] , {x,0,15}] 

and a plot similar to that of Figure 1.1 will be produced. Similarly, we can 

plot J n (x), for n = 1,2,3 Note that the graph of Jq{x) suggests that it has 

infinitely many positive zeros. 

On the open interval < x < oo, Bessel's equation has a two-dimensional 
space of solutions. However, it turns out that when p is a nonnegative integer, a 
second solution, linearly independent from the Bessel function of the first kind, 

3 For a very brief introduction to Mathematica, the reader can refer to Appendix A. 
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Figure 1.2: Graph of the Bessel function J\{x). 



cannot be obtained directly by the generalized power series method that we have 
presented. To obtain a basis for the space of solutions, we can, however, apply 
the method of variation of parameters just as we did in the previous section for 
the Cauchy-Euler equation; namely, we can set 

y = v(x)J p (x), 

substitute into Bessel's equation and solve for v(x). If we were to carry this out 
in detail, we would obtain a second solution linearly independent from J p (x). 
Appropriately normalized, his solution is often denoted by Y p (x) and called the 
p-th Bessel function of the second kind. Unlike the Bessel function of the first 
kind, this solution is not well-behaved near x = 0. 

To see why, suppose that yi(x) and 2/2(2:) is a basis for the solutions on the 
interval < x < 00, and let 1^(2/1,2/2) be their Wronskian, defined by 



W(y 1 ,y 2 ){x) 



2/1 (x) y[(x) 
2/2 (x) 2/ 2 (x) 



This Wronskian must satisfy the first order equation 

d 



dx 



(xW( yi ,y 2 )(x))=0, 



as one verifies by a direct calculation: 

£^(x2/i2/ 2 - xy 2 y[) = yix-^-(xy' 2 ) - y 2 x^-(xy' 1 ) 



= -(x 2 - n 2 )(y 1 y 2 - 2/22/1) = 0. 



2G 



Thus 

c 

xW(yi,y 2 ){x) = c, or W{y 1 ,y 2 ){x) = -, 

where c is a nonzero constant, an expression which is unbounded as x — > 0. 
It follows that two linearly independent solutions yi(x) and y 2 (x) to Bessel's 
equation cannot both be well-behaved at x = 0. 

Let us summarize what we know about the space of solutions to Bessel's 
equation in the case where p is an integer: 

• There is a one-dimensional space of real analytic solutions to (1.17), which 
are well-behaved as x — > 0. 

• This one-dimensional space is generated by a function J p (x) which is given 
by the explicit power series formula 

m=0 v ^ ; 



Exercises: 

1.4.1. Using the explicit power series formulae for Jo(x) and Ji(x) show that 

— Jo(x) = —J\(x) and —(xJ\(x)) — xJo(x). 
ax ax 

1.4.2. The differential equation 

is sometimes called a modified Bessel equation. Find a generalized power series 
solution to this equation in the case where p is an integer. (Hint: The power 
scries you obtain should be very similar to the power series for J p (x).) 

1.4.3. Show that the functions 

yi(x) = — \= cos a; and y 2 (x) = -^= sin x 
\Jx ^Jx 

are solutions to Bessel's equation 



*4 (x^)+(x 2 -p 2 )y = 0, 



dx \ dx 

in the case where p = 1/2. Hence the general solution to Bessel's equation in 
this case is 

1 1 

y = C\ —j= cos x + c 2 —= sin x. 



27 



1.4.4. To obtain a nice expression for the generalized power series solution to 
Bessel's equation in the case where p is not an integer, it is convenient to use 
the gamma function defined by 




- x e- l dt. 



a. Use integration by parts to show that T(x + 1) = xT(x). 

b. Show that T(l) = 1. 

c. Show that 

r(n + l) = n\ = n{n- l)---2- 1, 
when n is a positive integer. 



and use the recursion formula (1.21) to obtain the following generalized power 
series solution to Bessel's equation (1.17) for general choice of p: 



d. Set 



1 



a = 



2pr(p + i) 



y = W = (f) P E(- 1 ) 



m 



m!T(p + m+ 1) V2 




m=0 
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Chapter 2 



Symmetry and 
Orthogonality 



2.1 Eigenvalues of symmetric matrices 

Before proceeding further, we need to review and extend some notions from 
vectors and matrices (linear algebra), which the student should have studied 
in an earlier course. In particular, we will need the amazing fact that the 
eigenvalue-eigenvector problem for an n x n matrix A simplifies considerably 
when the matrix is symmetric. 

An n x n matrix A is said to be symmetric if it is equal to its transpose A T . 
Examples of symmetric matrices include 



Alternatively, we could say that an n x n matrix A is symmetric if and only if 



for every choice of n-vectors x and y. Indeed, since x • y = x T y, equation (2.1) 
can be rewritten in the form 



which holds for all x and y if and only if A = A T . 

On the other hand, an n x n real matrix B is orthogonal if its transpose is 
equal to its inverse, B T = B^ 1 . Alternatively, an n x n matrix 




and 




x • (Ay) = (Ax) ■ y. 



(2.1) 



x T Ay = (Ax) T y = x T A T y, 



B = (bib 2 "-b n ) 
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is orthogonal if its column vectors bi, b 2 , ... , b„ satisfy the relations 

bi-bi = l, bi-b 2 = 0, bi-b„ = 0, 

b 2 • b 2 = 1, • • • , b 2 • b„ = 0, 

b„ • b„ = 1. 

Using this latter criterion, we can easily check that, for example, the matrices 
/ a ■ q\ / 1/3 2/3 2/3N 

(=:-/)■ - {-iir-tt, 

are orthogonal. Note that since 

B T B = 1 (det S) 2 = (det S T )(dct B) = dct(B T B) = 1, 

the determinant of an orthogonal matrix is always ±1. 

Recall that the eigenvalues of an n x n matrix A are the roots of the poly- 
nomial equation 

det (4 - XI) = 0. 

For each eigenvalue Aj, the corresponding eigenvectors are the nonzero solutions 
x to the linear system 

(A - A7)x = 0. 

For a general n x n matrix with real entries, the problem of finding eigenvalues 
and eigenvectors can be complicated, because eigenvalues need not be real (but 
can occur in complex conjugate pairs) and in the "repeated root" case, there 
may not be enough eigenvectors to construct a basis for R n . We will see that 
these complications do not occur for symmetric matrices. 

Spectral Theorem. 1 Suppose that A is a symmetric n x n matrix with real 
entries. Then its eigenvalues are real and eigenvectors corresponding to distinct 
eigenvectors are orthogonal. Moreover, there is an n x n orthogonal matrix B 
of determinant one such that B~ x AB = B T AB is diagonal. 

Sketch of proof: The reader may want to skip our sketch of the proof at first, 
returning after studying some of the examples presented later in this section. 
We will assume the following two facts, which are proven rigorously in advanced 
courses on mathematical analysis: 

1. Any continuous function on a sphere (of arbitrary dimension) assumes its 
maximum and minimum values. 

2. The points at which the maximum and minimum values are assumed can 
be found by the method of Lagrange multipliers (a method usually dis- 
cussed in vector calculus courses). 



1 This is called the "spectral theorem" because the spectrum is another name for the set of 
eigenvalues of a matrix. 
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The equation of the sphere 5™ 1 in R" is 

x i + x \ + ' ' ' + x n = 1; or x T x = 1, where x = 



( Xl \ 

\X„J 



We let 

g(x) = g(x ll x 2l ... , x n ) = x T x - 1, 

so that the equation of the sphere is given by the constraint equation <?(x) = 0. 
Our approach consists of finding of finding the point on S 1 " -1 at which the 
function 

/( x ) = f(xi,x 2 , ... ,x n ) = x T Ax 

assumes its maximum values. 

To find this maximum using Lagrange multipliers, we look for "critical 
points" for the function 

H(x,X) = H(x!,x 2 , ... , x n , A) = /(x) - A.g(x). 

These are points at which 

V/(a;i,x 2 , . . . ,x n ) = Wg(xi,x 2 , ■ ■ ■ ,x n ), and g(xi, x 2 , ■ ■ ■ , x n ) = 0. 

In other words, these are the points on the sphere at which the gradient of / is 
a multiple of the gradient of g, or the points on the sphere at which the gradient 
of / is perpendicular to the sphere. 
If we set 

(an ai2 • »in\ 
a 2 i a 22 ■ a 2r , 

\a n i o, n2 ■ a nn J 
a short calculation shows that at the point where / assumes its maximum, 

dxi 

or equivalcntly, 



2a il x 1 + 2a l2 x 2 



2\ Xi = 0, 



Ax - Ax = 0. 



We also obtain the condition 



dH 



-S(x) = 0, 



which is just our constraint. Thus the point on the sphere at which / assumes 
its maximum is a unit-length eigenvector t>i, the eigenvalue being the value Ai 
of the variable A. 

Let W be the "linear subspace" of R n defined by the homogeneous linear 
equation t>i • x = 0. The intersection S""" 1 n W is a sphere of dimension n — 2. 
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We next use the method of Lagrange multipliers to find a point on S™ 1 n W 
at which / assumes its maximum. To do this, we let 

5i(x)=x T x-l, g 2 (x) = bi-x. 

The maximum value of / on S" 1-1 n W will be assumed at a critical point for 
the function 

H(x, X,fi) = /(x) - A. 9l (x) - fig 2 (x). 
This time, differentiation shows that 

— = 2anxi + 2a i2 x 2 H h 2a in x n - 2Xx { - \ibi = 0, 

or equivalently, 

Ax - Ax - jtibi = 0. (2.2) 
It follows from the constraint equation b x • x = that 
bi • (Ax) = bl{Ax) = (bjA)x = (bf A T )x 

= (Abi) T x = (Aibi) T x = Aibi • x = 0. 

Hence it follows from (2.2) that Ax - Ax = 0. Thus if b 2 is a point on S*™" 1 n W 
at which / assumes its maximum, b 2 must be a unit-length eigenvector for A 
which is perpendicular to bi. 

Continuing in this way we finally obtain n mutually orthogonal unit-length 
eigenvectors b l5 b 2 , . . . , b„. These eigenvectors satisfy the equations 

Abi = Aibi, ,4b2 = A2b2, ... Ab n = A 2 b„, 
which can be put together in a single matrix equation 

(Abi Ab 2 ■ Ab n ) = (Aibi A 2 b 2 • A„b„) , 
or equivalently, 












A (bi b 2 • b„) = (bi b 2 • b„) 





A 2 















If we set 

B = (bi b 2 • b„) , 

this last equation becomes 












AB = B 





A 2 
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Of course, B is an orthogonal matrix, so it is invertible and we can solve for A, 
obtaining 



A = B 







A 2 





B-\ 


or B~ X AB = 


A 





A 2 







\° 











u 





A J 



We can arrange that the determinant of -B be one by changing the sign of one 
of the columns if necessary. 

A more complete proof of the theorem is presented in more advanced courses 
in linear algebra. 2 In any case, the method for finding the orthogonal matrix B 
such that B T AB is diagonal is relatively simple, at least when the eigenvalues 
are distinct. Simply let B be the matrix whose columns are unit-length eigen- 
vectors for A. In the case of repeated roots, we must be careful to choose a 
basis of unit-length eigenvectors for each eigenspace which are perpendicular to 
each other. 

Example. The matrix 




is symmetric, so its eigenvalues must be real. Its characteristic equation is 



5 - A 4 
4 5 - A 
1 - A 



[(A-5) 2 -16](A-l) 



= (A 2 - 10A + 9)(A - 1) = (A - 1) 2 (A - 9), 

which has the roots Ai — 1 with multiplicity two and A 2 — 9 with multiplicity 
one. 

Thus we are in the notorious "repeated root" case, which might be expected 
to cause problems if A were not symmetric. However, since A is symmetric, 
the Spectral Theorem guarantees that we can find a basis for J? 3 consisting of 
eigenvectors for A even when the roots are repeated. 

We first consider the eigenspace W\ corresponding to the eigenvalue Ai = 1, 
which consists of the solutions to the linear system 



or 



5-1 4 







h 


4 5-1 





)( 


b 2 





1 - 1 




b 3 


4&i 


+46 2 


= 0, 




4&i 


+46 2 


= o, 









= 0. 





2 There are many excellent linear algebra texts that prove this theorem in detail; one good 
reference is Bill Jacob, Linear algebra, W. H. Freeman, New York, 1990; see Chapter 5. 
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The coefficient matrix of this linear system is 




Applying the elementary row operations to this matrix yields 



4 


4 






1 — 1 




{ 1 1 





4 


4 




-I 


4 




+ 00 











o / 







o / 


loo 






Thus the linear system is equivalent to 

&i + & 2 = 0, 
=0, 
=0. 

Thus W\ is a plane with equation b\ + b 2 = 0. 

We need to extract two unit length eigenvectors from W\ which are perpen- 
dicular to each other. Note that since the equation for W\ is b\ + 62 = 0, the 
unit length vector 

1 

± 


is perpendicular to W\. Since 

b x = I I 6W1, we find that b 2 = n x h x = \ ^ | e Wi. 

The vectors bi and b 2 arc unit length elements of W\ which are perpendicular 
to each other. 

Next, we consider the eigenspace Wg corresponding to the eigenvalue A2 = 9, 
which consists of the solutions to the linear system 





or 

-4&i +46 2 = 0, 

46i -46 2 = 0, 



-003 

The coefficient matrix of this linear system is 



0. 



-4 


4 





4 


-4 











—I 
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Applying the elementary row operations to this matrix yields 



-4 


4 





4 


-4 











-i 




Thus the linear system is equivalent to 

h - b 2 

h 



and we sec that 



0, 
0, 

o, 



Wo 



span 



We set 



b 3 = 




Theory guarantees that the matrix 
B = 




whose columns arc the eigenvectors h»i, h»i, and b 3 , will satisfy 

B~ 1 AB = 




Moreover, since the eigenvectors we have chosen are of unit length and perpen- 
dicular to each other, the matrix B will be orthogonal. 

Exercises: 

2.1.1. Find a 2 x 2-matrix B such that B is orthogonal and B~ 1 AB is diagonal, 
where 

' 5 4 
4 5 

2.1.2. Find a 2 x 2-matrix B such that B is orthogonal and B~ x AB is diagonal, 
where 

' 3 4 
4 9 
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2.1.3. Find a 3 x 3-matrix B such that B is orthogonal and B 1 AB is diagonal, 
where 

/ 5 4 \ 
A= 4 5 . 
V 1 j 

2.1.4. Find a 3 x 3-matrix B such that B is orthogonal and B~ x AB is diagonal, 
where 

/ — 1 2 \ 
A= 2 2. 
V 2 1 j 

2.1.5. Show that the nxn identity matrix / is orthogonal. Show that if B\ and 
B 2 are nxn orthogonal matrices, so is the product matrix BiB 2 , as well as the 
inverse B7 1 . 



2.2 Conic sections and quadric surfaces 

The theory presented in the previous section can be used to "rotate coordinates" 
so that conic sections or quadric surfaces can be put into "canonical form." 

A conic section is a curve in R 2 which is described by a quadratic equation, 
such as the equation 

ax\ + 2bx\X2 + cx\ = 1, (2.3) 
where a, b and c are constants, which can be written in matrix form as 



If we let 




we can rewrite (2.3) as 

x T ^x=l, (2.4) 

where A is a symmetric matrix. 

According to the Spectral Theorem from Section 2.1, there exists a 2 x 2 
orthogonal matrix B of determinant one such that 

B^AB = B T AB = 

If we make a linear change of variables, 

x = Sy, 



/ Ax \ 
I A 2 ■ 
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) 




o' 


A 2 J 


S 



then since x T = y T B T , equation (2.3) is transformed into 

y T {B T AB)y =1, ( Wl y 2 
or equivalently, 

\ iy l + \ 2V l = l. (2.5) 
In the new coordinate system (2/1,2/2), it is easy to recognize the conic section: 

• If Ai and A2 are both positive, the conic is an ellipse. 

• If Ai and A 2 have opposite signs, the conic is an hyperbola. 

• If Ai and A2 are both negative, the conic degenerates to the the empty set, 
because the equation has no real solutions. 

In the case where Ai and A 2 are both positive, we can rewrite (2.5) as 



(v/i/at) 2 (yiA^ 

from which we recognize that the semi-major and semi-minor axes of the ellipse 
are y/l/Xi and ^l/Aa. 

The matrix B which relates x and y represents a rotation of the plane. To 
see this, note that the first column b>i of B is a unit-length vector, and can 
therefore be written in the form 

. / cos 9 

for some choice of 9. The second column b 2 is a unit-vector perpendicular to 
bi and hence 

- sin V s 



cos 9 

We must take the plus sign in this last expression, because det B = 1 . Thus 

" G) - (s# ■ * (!) - (■«/ 

or equivalently, B takes the standard basis vectors for R 2 to vectors which have 
been rotated counterclockwise through an angle 9. By linearity, 



B = 



cos 9 — sin 9 
sin 9 cos 9 



must rotate every element of J? 2 counterclockwise through an angle 9. Thus once 
we have sketched the conic in (j/i, ^-coordinates, we can obtain the sketch in 
(a;i, a^-coordinates by simply rotating counterclockwise through the angle 9. 
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Example. Let's consider the conic 

hx\ — 6x1X2 + hx\ — 1, (2.6) 

or equivalently, 

The characteristic equation of the matrix 
is 

(5 — A) 2 — 9 = 0, A 2 -10A + 16 = 0, (A - 2) (A - 8) = 0, 

and the eigenvalues are Ai = 2 and A 2 = 8. Unit-length eigenvectors corre- 
sponding to these eigenvalues are 



bi 



1/V2 \ ( -1/V2 



The proof of the theorem of Section 2.1 shows that these vectors are orthogonal 
to each other, and hence the matrix 



1/V2 -1/V2 



B = 



is an orthogonal matrix such that 

Note that B represents a counterclockwise rotation through 45 degrees. If 
we define new coordinates (2/1,2/2) by 

xi \ j; I IJl 



x 2 I V 2/2 



equation (2.6) will simplify to 



(»• *)(2!!)(S)^ 

or 

2 yl + gy 2 2 = — + — *f=- = 1. 

1 2 (1/V2) 2 (i/Vs) 2 

We recognize that this is the equation of an ellipse. The lengths of the semimajor 
and scmiminor axes are l/y/2 and 1/(2^2). 



; J ,8 





Figure 2.1: Sketch of the conic section 5x\ ~ 6x1X2 + 5x 2 — 1 = 0. 

The same techniques can be used to sketch quadric surfaces in R 3 , surfaces 
defined by an equation of the form 

(an a 12 a 13 
a.21 a 2 2 a 23 
^31 «32 033 

where the ay's are constants. If we let 

(an a\2 ai3 
021 a 2 2 a 2 3 
^31 a 32 «33 

we can write this in matrix form 

x T Ax=l. (2.7) 

According to the Spectral Theorem, there is a 3 x 3 orthogonal matrix B of 
determinant one such that B T AB is diagonal. We introduce new coordinates 

y = 2/2 by setting x = By, 

w 

and equation (2.7) becomes 

y T (B T AB)y = 1. 

Thus after a suitable linear change of coordinates, the equation (2.7) can be put 
in the form 

■ Aj \ / yi 

(2/12/22/3)! A 2 2/2 ) = 1, 

A 3 I \ y 3 
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Figure 2.2: An ellipsoid. 



or 

MVi + A22/2 + A 32/3 = !) 
where Ai, A2, and A3 are the eigenvalues of A. It is relatively easy to sketch the 
quadric surface in the coordinates (2/1,2/2,2/3)- 

If the eigenvalues are all nonzero, we find that there are four cases: 

• If Ai, A2, and A3 are all positive, the quadric surface is an ellipsoid. 

• If two of the A's are positive and one is negative, the quadric surface is an 
hyperboloid of one sheet. 

• If two of the A's are negative and one is positive, the quadric surface is an 
hyperboloid of two sheets. 

• If Ai, A2, and A3 are all negative, the equation represents the empty set. 

Just as in the case of conic sections, the orthogonal matrix B of determinant 
one which relates x and y represents a rotation. To see this, note first that since 
B is orthogonal, 

(Bx) ■ (By) = x T S T By = x T /y = x ■ y. (2.8) 

In other words, multiplication by B preserves dot products. It follows from this 
that the only real eigenvalues of B can be ±1. Indeed, if x is an eigenvector for 
B corresponding to the eigenvalue A, then 

A 2 (x • x) = (Ax) • (Ax) = (Bx) • (Bx) = x • x, 
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so division by x • x yields A 2 = 1 . 

Since detB = 1 and det-B is the product of the eigenvalues, if all of the 
eigenvalues are real, A = 1 must occur as an eigenvalue. On the other hand, 
non-real eigenvalues must occur in complex conjugate pairs, so if there is a non- 
real eigenvalue fi + iv, then there must be another non-real eigenvalue /i — iv 
together with a real eigenvalue A. In this case, 

1 = det B = AO + iv){fi - iv) = A(> 2 + v 2 ). 

Since A = ±1, we conclude that A = 1 must occur as an eigenvalue also in this 
case. 

In either case, A — 1 is an eigenvalue and 

W x = {x g R 3 : Bx = x} 

is nonzero. It is easily verified that if dimWi is larger than one, then B must 
be the identity. Thus if B ^ 7, there is a one-dimensional subspace W\ of R 3 
which is left fixed under multiplication by B. This is the axis of rotation. 

Let Wi denote the orthogonal complement to W\ . If x is a nonzero element 
of Wi and y g W^, it follows from (2.8) that 

(By) • x = (By) ■ (Bx) = y T B T Bx = y T /x = y • x = 0, 

so By g Wi~. Let y, z be elements of such that 

y -y = 1, yz = 0, z • z = 1; 

we could say that {y,z} form an orthonormal basis for W±. By (2.8), 

(By) ■ (By) = 1, (By) ■ (Bz) = 0, (Bz) ■ (Bz) = 1. 

Thus By must be a unit-length vector in and there must exist a real number 
9 such that 

By = cos Oy + sin Bz. 

Moreover, Bz must be a unit-length vector in which is perpendicular to 
By and hence 

Bz = ±(— sin#y + cos#z). 
However, we cannot have 

Bz = — (— sin#y + cos#z). 

Indeed, this would imply (via a short calculation) that the vector 

u = cos(fl/2)y + sin(0/2)z 

is fixed by B, in other words Bu = u, contradicting the fact that u g . 
Thus we must have 

Bz — — sin 9y + cos 6z, 
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Figure 2.3: Hyperboloid of one sheet. 



and multiplication by B must be a rotation in the plane through an angle 
9. Moreover, it is easily checked that y + iz and y — iz are eigenvectors for B 
with eigenvalues 



We can therefore conclude that a 3 x 3 orthogonal matrix B of determinant 
one represents a rotation about an axis (which is the cigenspacc for eigenvalue 
one) and through an angle 9 (which can be determined from the eigenvalues of 
B, which must be 1 and e ±lB ). 

Exercises: 

2.2.1. Suppose that 



= cos 9 ± i sin 9. 




a. Find an orthogonal matrix B such that B T AB is diagonal. 

b. Sketch the conic section 2xf + 6x1X2 + 2x\ = 1. 

c. Sketch the conic section 2(xi — l) 2 + 6(2; 1 — l)x% + 2x\ = 1. 



2.2.2. Suppose that 
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2 

Figure 2.4: Hyperboloid of two sheets. 



a. Find an orthogonal matrix B such that B T AB is diagonal. 

b. Sketch the conic section 2x\ — 4a; 1X2 + 5x\ = 1. 

c. Sketch the conic section 2xf — Ax\X2 + 5x\ — Ax\ + 4x2 = - 



1. 



2.2.3. Suppose that 



A = 



( 



4 2 
2 1 



) 



a. Find an orthogonal matrix B such that B T AB is diagonal. 



b. Sketch the conic section 4xf + 4x\X2 + x\ — \fhx\ + 2\/hx2 = 0. 

2.2.4. Determine which of the following conic sections are ellipses, which are 
hyperbolas, etc.: 

a. x\ + Ax\X2 + 3^2 = 1. 

b. x\ + 6x1X2 + lOxl = 1. 

c. — 3xf + 6x1X2 — 4^2 = 1. 

d. — x\ + Axix 2 — ix\ = 1. 

2.2.5. Find the semi-major and semi-minor axes of the ellipse 



5x1 + 6x1X2 + 5x2 



= 4. 
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2.2.6. Suppose that 

/ 2 \ 
A= 2 3 . 
\ 9 / 

a. Find an orthogonal matrix B such that B T AB is diagonal. 

b. Sketch the quadric surface Ax\x 2 + 3x 2 + 9xg = 1. 

2.2.7. Determine which of the following quadric surfaces are ellipsoids, which 
are hyperboloids of one sheet, which are hyperboloids of two sheets, etc.: 

a. x\ — x\ — x\ — 6x2X3 = 1. 

b. x\ + x\ — x\ — 6x1X2 = 1. 

c. x\ + x\ + x\ + 4xix 2 + 2x 3 = 0. 

2.2.8. a. Show that the matrix 

/ 1/3 2/3 2/3\ 
B = -2/3 -1/3 2/3 
V 2/3 -2/3 1/3/ 

is an orthogonal matrix of dctcminant one. Thus multiplication by B is rotation 
about some axis through some angle. 

b. Find a nonzero vector which spans the axis of rotation. 

c. Determine the angle of rotation. 



2.3 Orthonormal bases 

Recall that if v is an element of K 3 , we can express it with respect to the 
standard basis as 

v = ai + bj + ck, where a = v • i, b = v • j, c = v-k. 

We would like to extend this formula to the case where the standard basis 
{i,j,k} for R 3 is replaced by a more general "orthonormal basis." 

Definition. A collection of n vectors bi, b 2 , ... , b„ in R n is an orthonormal 
basis for R n if 



bi-bi = l, bi-b 2 = 0, bi-b„ = 0, 

b 2 • b 2 = 1, • • • , b 2 • b„ = 0, 



(2.9) 



b„ • b„ = 1. 
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From the discussion in Section 2.1, we recognize that the term orthonormal 
basis is just another name for a collection of n vectors which form the columns 
of an orthogonal n x n matrix. 

It is relatively easy to express an arbitrary vector f € R n in terms of an 
orthonormal basis bi, b 2 , . . . , b„: to find constants ci, c 2 , . . . , c n so that 

f = cibi + c 2 b 2 H h (2.10) 

we simply dot both sides with the vector bj and use (2.9) to conclude that 

Ci = bj • f . 

In other words, if bi, b 2 , . . . , b„ is an orthonormal basis for R n , then 

fen" =* f = (f ■bi)bi + --- + (f-b n )b n) (2.11) 

a generalization of the formula we gave for expressing a vector in R 3 in terms 
of the standard basis {i,j,k}. 

This formula can be helpful in solving the initial value problem 



dx 

~dt 



= Ax, 



x(0) = f , 



(2.12) 



in the case where A is a symmetric n x n matrix and f is a constant vector. 
Since A is symmetric, the Spectral Theorem of Section 2.1 guarantees that the 
eigenvalues of A are real and that we can find an n x n orthogonal matrix B 
such that 





fXi 


• 




B~ 1 AB = 





A 2 • 


• 




V° 


• 


• KJ 



If we set x — By, then 



B-f- = — = Ax = ABy 
dt dt 



dy 

dt 



= B- x ABy. 



Thus in terms of the new variable 



2/2 

\VnJ 



we have 



fdyi/dt\ 
dyijdt 

\dy n /dtj 



Ai 

A 2 
\o 



\ 






2/2 


/ 


w 



so that the matrix differential equation decouples into n noninteracting scalar 
differential equations 

dyi/dt = Aiyi, 
dyijdt = X 2 y2, 



dy n /dt = X n y n . 
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We can solve these equations individually, obtaining the general solution 



(Vi\ 




/ Cl e Al *\ 


V2 




c 2 e A2 * 


\Vn) 







y = 



where ci, c 2 , ... , c„ are constants of integration. Transferring back to our 
original variable x yields 



x = B 



/ Cl e Al *\ 
c 2 e A2 * 



= cib ie Al * + c 2 b 2 e A2t + • • • + c„b„e 



A 2 t 



,A„t 



\c„e A "7 

where bi, b 2 , . . . , b„ are the columns of B. Note that 

x(0) = cibi + c 2 b 2 H h c„b„. 

To finish solving the initial value problem (2.12), we need to determine the 
constants c\, c 2 , . . . , c„ so that 

cibi + c 2 b 2 + h c„b„ = f . 

It is here that our formula (2.11) comes in handy; using it we obtain 

x = (f • b!)b ie Al * + • • • + (f • b„)b n e A " 4 . 

Example. Suppose we want to solve the initial value problem 



~dt 



<h 4 CP 
= I 4 5 | x. 

.0 1; 



x(0) = f , where 



We saw in Section 2.1 that A has the eigenvalues Ai — 1 with multiplicity two 
and A 2 = 9 with multiplicity one. Moreover, the orthogonal matrix 




has the property that 



B 



1 
B~ 1 AB= I 1 



9 

Thus the general solution to the matrix differential equation dx./dt — Ax is 



cie" 

x = B\ C2 e 4 



\C3e 
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cibie* + c 2 b 2 e* + c 3 b 3 e 9t , 
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where 



bi 





Setting t = in our expression for x yields 

x(0) = cibi + c 2 b 2 + c 3 b 3 . 
To solve the initial value problem, we employ (2.11) to see that 



Cl 




C2 




V2, 



c 3 




2V2. 



Hence the solution is 




o9* 



Exercise: 

2. 3.1. a. Find the eigenvalues of the symmetric matrix 

/5 4 0\ 
. _ 4 5 

4 2' 
\0 2 I J 

b. Find an orthogonal matrix B such that B~ x AB is diagonal. 

c. Find an orthonormal basis for R consisting of eigenvectors of A. 

d. Find the general solution to the matrix differential equation 

dx 

— = Ax. 
dt 

e. Find the solution to the initial value problem 

dx 



dt 



= Ax, 



x(0) 



3 


V 2 / 
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2. 3. 2. a. Find the eigenvalues of the symmetric matrix 



A = 



-3 



2 

-4 
2 







J 



(Hint: To find roots of the cubic, try A = 1 and A = — 1.) 

b. Find an orthogonal matrix B such that B~ 1 AB is diagonal. 

c. Find an orthonormal basis for J? 3 consisting of eigenvectors of A. 



d. Find the general solution to the matrix differential equation 



— = Ax. 
dt 



e. Find the solution to the initial value problem 




dx 



x(0) = 2 



:) 



2.3.3. (For students with access to Mathematica) a. Find the eigenvalues of the 
matrix 



by running the Mathematica program 

a = {{2, 1,1}, {1,3,2}, {1,2, 4}}; Eigenvalues [a] 

The answer is so complicated because Mathematica uses exact but complicated 
formulae for solving cubics discovered by Cardan and Tartaglia, two sixteenth 
century Italian mathematicians. 

b. Find numerical approximations to these eigenvalues by running the program 
Eigenvalues [N [a] ] 

The numerical values of the eigenvalues are far easier to use. 

c. Use Mathematica to find numerical values for the eigenvectors for A by 
running the Mathematica program 

Eigenvectors [N [a] ] 

and write down the general solution to the matrix differential equation 
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Figure 2.5: Two carts connected by springs and moving along a friction- free 
track. 



2.4 Mechanical systems 

Mechanical systems consisting of weights and springs connected together in an 
array often lead to initial value problems of the type 

x(0) = f, f = 0, (2.13) 

where A is a symmetric matrix and f is a constant vector. These can be solved 
by a technique similar to that used in the previous section. 

For example, let us consider the mechanical system illustrated in Figure 2.5. 
Here we have two carts moving along a friction-free track, each containing mass 
m and attached together by three springs, with spring constants ki, k 2 and k%. 
Let 

x\{t) = the position of the first cart to the right of equilibrium, 
x 2 (t) = the position of the second cart to the right of equilibrium, 
Fi = force acting on the first cart, 
F 2 = force acting on the second cart, 

with positive values for F\ or F 2 indicating that the forces are pulling to the 
right, negative values that the forces are pulling to the left. 

Suppose that when the carts are in equilibrium, the springs are also in equi- 
librium and exert no forces on the carts. In this simple case, it is possible to 
reason directly from Hooke's law that the forces F\ and F 2 must be given by 
the formulae 

Fi = -k\Xi + k 2 (x 2 - x\), F 2 = k 2 (xi - x 2 ) - k 3 x 2 , 

but it becomes difficult to determine the forces for more complicated mechanical 
systems consisting of many weights and springs. Fortunately, there are some 
simple principles from physics which simply the procedure for finding the forces 
acting in such mechanical systems. 
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The easiest calculation of the forces is based upon the notion of work. On 
the one hand, the work required to pull a weight to a new position is equal to 
the increase in potential energy imparted to the weight. On the other hand, we 
have the equation 

Work = Force x Displacement, 

which implies that 

Work Change in potential energy 



Force = 



Displacement Displacement 



Thus if V{x\, X2) is the potential energy of the configuration when the first cart 
is located at the point x\ and the second cart is located at the point x 2 , then 
the forces are given by the formulae 

F 1 = -™ F 2 = -™ 

± 1 o ' 2 n 

OX\ OX2 

In our case, the potential energy V is the sum of the potential energies stored 
in each of the three springs, 

V(x 1 ,x 2 ) = \^k\x\ + ^k 2 {x 1 - x 2 ) 2 + ^k 3 xl, 

and hence we obtain the formulae claimed before: 

dV 

F 1 = -—— = -k 1 x 1 + k 2 (x 2 - xi), 

OX\ 

dV 

F 2 = - ^ — = k 2 (xi - x 2 ) - k 3 x 2 . 
ox 2 

It now follows from Newton's second law of motion that 
Force = Mass x Acceleration, 

and hence 

d 2 X! d 2 x 2 
Thus we obtain a second-order system of differential equations, 

d 2 X! 

m ~^2 = -fci^i + k 2 (x 2 - x x ) = -(fci + k 2 )x\ + k 2 x 2 , 
d 2 x 2 

m ~Tj-r = k 2 {xi - x 2 ) - k 3 x 2 = k 2 xi - (k 2 + k 3 )x 2 . 
We can write this system in matrix form as 
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Note that A is indeed a symmetric matrix. The potential energy is given by the 
expression 



V(x 1 ,x 2 ) = -\ (xix 2 ) A r 



Example. Let us consider the special case of the mass-spring system in which 

m = k\ = k 2 = k% = 1, 
so that the system (2.14) becomes 

To find the eigenvalues, we must solve the characteristic equation 
det(~ 2 ~ A I _ A ) =(A + 2) 2 -l = 0, 

which yields 

A = -2 ± 1. 

The eigenvalues in this case are 

Ai = — 1, and A 2 = —3. 
The eigenspace corresponding to the eigenvalue —1 is 

W-i = {b e R 2 : (A + I)b = 0} = . . . = span ^ j 

It follows from the argument in Section 2.1 that the eigenspace corresponding 
to the other eigenvalue is just the orthogonal complement 



W-3 = span 



-1 
1 



Unit length eigenvectors lying in the two eigenspaces are 



h _ [ l/\/2 \ , _ ( -l/x/2 



The theorem of Section 2.1 guarantees that the matrix 

B = 



l/s/2 -1/V2 
l/y/2 l/y/2 J' 

whose columns are hi and b 2 , will diagonalizc our system of differential equa- 
tions. 



51 



Indeed, if we define new coordinates (2/1,2/2) by setting 

xi\ = ( 1/V2 -1/V2 \( yi 
x 2 ) \ 1/^/2 1/V2 ) \ y 2 

our system of differential equations transforms to 

d 2 y 1 /dt 2 = -yi, 
d 2 y 2 /dt 2 = -32/ 2 . 

We set uji = 1 and uj 2 = \/3, so that this system assumes the familiar form 

d 2 yi /dt 2 + iv 2 Vl = 0, 
d 2 y 2 /dt 2 + uo 2 y 2 = 0, 

a system of two noninteracting harmonic oscillators. 
The general solution to the transformed system is 

2/1 = ai cosujit + b\ sincjii, y 2 — a 2 cos u> 2 t + b 2 sintv 2 t. 

In the original coordinates, the general solution to (2.15) is 

x\ \ _ / l/\/2 —1/V2 \ / ai coswit + 61 sin 
x 2 / \ l/\/2 l/\/2 / \ a 2 cos u 2 t + b 2 sin uj 2 t 

or equivalently, 

x = h>i(ai coscjit + bi sinu>ii) + b 2 (a2 cos^i + b 2 sina> 2 i)- 

The motion of the carts can be described as a general superposition of two 
modes of oscillation, of frequencies 

U>! LU 2 

— and — . 
2tt 2tt 

Exercises: 

2. 4.1. a. Consider the mass-spring system with two carts illustrated in Figure 2.5 
in the case where fci = 4 and m = k 2 = k 3 = 1. Write down a system of second- 
order differential equations which describes the motion of this system. 

b. Find the general solution to this system. 

c. What are the frequencies of vibration of this mass-spring system? 

2. 4. 2. a. Consider the mass-spring system with three carts illustrated in Fig- 
ure 2.5 in the case where m = k\ = k 2 = = fej = 1. Show that the motion of 
this system is described by the matrix differential equation 

d 2 * 

di 2 









1 





Ax, 


where 




-2 








\ 


1 
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Figure 2.6: Three carts connected by springs and moving along a friction-free 
track. 



b. Find the eigenvalues of the symmetric matrix A. 

c. Find an orthonormal basis for R 3 consisting of eigenvectors of A. 

d. Find an orthogonal matrix B such that B~ 1 AB is diagonal. 

e. Find the general solution to the matrix differential equation 

rf 2 x 



dt 2 



= Ax. 



f. Find the solution to the initial value problem 
d 2 x 



dt 2 



Ax, 



x(0) 



~dt 



(0) 



2. 4. 3. a. Find the eigenvalues of the symmetric matrix 



A = 



b. What are the frequencies of oscillation of a mechanical system which is 
governed by the matrix differential equation 



r 2 


1 





M 


i 


-2 


1 








1 


-2 


1 


\o 





1 


-v 



dt 2 



= Ax? 



2.5 Mechanical systems with many degrees of 
freedom* 

Using an approach similar to the one used in the preceding section, we can 
consider more complicated systems consisting of many masses and springs. For 
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Figure 2.7: A circular array of carts and springs. 



example, we could consider the box spring underlying the mattress in a bed. 
Although such a box spring contains hundreds of individual springs, and hence 
the matrix A in the corresponding dynamical system constains hundreds of rows 
and columns, it is still possible to use symmetries in the box spring to simplify 
the calculations, and make the problem of determining the "natural frequencies 
of vibration" of the mattress into a manageable problem. 

To illustrate how the symmetries of a problem can make it much easier to 
solve, let us consider a somewhat simpler problem, a system of n carts containing 
identical weights of mass m and connected by identical springs of spring constant 
k, moving along a circular friction- free track as sketched in Figure 2.6. 

We choose a positive direction along the track and let denote the dis- 
placement of the i-th cart out of equilibrium position in the positive direction, 
for 1 < i < n. The potential energy stored in the springs is 

V(xi,... ,x n ) = \k(x n ~xi) 2 +\k{x 2 -xi) 2 +\k(x3-x 2 ) 2 +- ■ .+\k{x n -x n ^i) 2 







/ 
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•• 
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We can write this as 



V(xi,... , x n ) = - ^ kx T Ax, 



where 





/ 


Xi 


\ 




( 


-2 


1 


•• 


1 


\ 






Xi 








1 


-2 


1 •• 


• 




X = 




%3 











1 


-2 ■■ 


• 






V 




J 




V 


1 





•• 


■ -2 


/ 



or equivalently as 



^ n n 

V(xi , . . . 3 Xn') — — ^ ^ ^ ^ ^ ^ O'ij^iXj i 
i=l 2=1 
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where denotes the (i, j)-component of the matrix A. 

Just as in the preceding section, the force acting on the i-th cart can be 
calculated as minus the derivative of the potential energy with respect to the 
position of the i-th cart, the other carts being held fixed. Thus for example, 



dV 



^ n ^ n n 



ljXj, 



the last step obtained by using the symmetry of A. In general, we obtain the 
result: 

dV 

dxi 



Fi 



k dijXj, 



3 = 1 



which could be rewritten in matrix form as 

F = kAx. 

On the other hand, by Newton's second law of motion, 

d 2 x 



(2.16) 



in 



dt 2 



Substitution into (2.16) yields 
d 2 x 



to—— = kAx or 
dt 2 



F. 



d 2 x k A 

172" = -Ax > 
dt 1 to 



(2.17) 



where A is the symmetric matrix given above. To find the frequencies of vi- 
bration of our mechanical system, we need to find the eigenvalues of the n x n 
matrix A. 

To simplify the calculation of the eigenvalues of this matrix, we make use of 
the fact that the carts are identically situated — if we relabel the carts, shifting 
them to the right by one, the dynamical system remains unchanged. Indeed, 
lew us define new coordinates (j/i, . . . , y n ) by setting 



xi = 2/2, 
or in matrix terms, 



X2 = 2/3, 



%n— 1 l)n. •) %n — Ul ■> 



x = Ty, where T - 



( 1 

1 





\ 1 



o \ 



1 

o / 



Then y satisfies exactly the same system of differential equations as x: 

d 2 y k 



dt 2 



TO 



-Ay. 



(2.18) 
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On the other hand, 

-Ax => T- - 



k 



ATy 



dt 2 m dt 2 m 

Comparison with (2.18) yields 

A = T^AT, or TA 



or 



AT. 



dt 2 



m 



-T~ l ATy. 



In other words the matrices A and T commute. 

Now it is quite easy to solve the eigenvector-eigenvalue problem for T . In- 
deed, if x is an eigenvector for T corresponding to the eigenvalue A, the compo- 
nents of x must satisfy the vector equation Tx = Ax. In terms of components, 
the vector equation becomes 



x 2 = Xxi, x 3 = Xx 2 , ■ . . ,x n = Aa; n _i, x x = \x n . 
Thus X3 = X 2 xi, x 4 = X 3 xi, and so forth, yielding finally the equation 

xi = \ n xi. 

Similarly, 

x 2 = X n x 2 , ... , x n = X n x n . 
Since at least one of the a^'s is nonzero, we must have 

A™ = 1. 



(2.19) 



(2.20) 



This equation is easily solved via Euler's formula: 



3 27ri/nyin _ 



1 = e 2 ™ (e 2 "/")" = 1, and similarly [(e 
for < j < n — 1. Thus the solutions to (2.20) are 

A = 77^ , for 0<j <n-l, where 77 = e 27Ti/n 



)T = i. 



(2.21) 



For each choice of j, we can try to find eigenvectors corresponding to rf . If 
we set x\ — 1, we can conclude from (2.19) that 



x 2 = if, X 3 = T] 



2j 



7) 



(n-l)j 



thereby obtaining a nonzero solution to the eigenvector equation. Thus for each 
j, < j < n — 1, we do indeed have an eigenvector 



1 

rf 

T] 



\ 



2j 



I 



for the eigenvalue if. Moreover, each eigenspace is one-dimensional. 
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In the dynamical system that we are considering, of course, we need to solve 
the eigenvalue-eigenvector problem for A, not for T. Fortunately, however, since 
A and T commute, the eigenvectors for T are also eigenvectors for A. Indeed, 
since AT = TA, 

T(Ae 3 ) = A(Tej) = Atfej) = rf(Aej), 

and this equation states that Aej is an eigenvector for T with the same eigen- 
value as Bj . Since the eigenspaces of T are all one-dimensional, Aej must be a 
multiple of e^; in other words, 

Aej = Xjej, for some number Xj. 

To find the eigenvalues Xj for A, we simply act on ej by A: we find that the 
first component of the vector Aej = Xjej is 

-2 + rf + r/™- 1 ^' = -2 + (e 2 ^'/" + e -**H/n) = _2 + 2 cos(27rj/n), 

where we have used Euler's formula once again. On the other hand, the first 
component of ej is 1, so we immediately conclude that 

Xj = -2 + 2cos(27rj/n). 

It follows from the familiar formula 

cos(2a) = cos 2 a — sin 2 a 

that 

Xj = -2 + 2[cos(7rj/n)] 2 - 2[sin(7rj/n)] 2 = -4[sin(7rj/n)] 2 . 

Note that Xj = X n -j, and hence the eigenspaces for A are two-dimensional, 
except for the eigenspace corresponding to j = 0, and if n is even, to the 
eigenspace corresponding to j = n/2. Thus in the special case where n is odd, 
all the eigenspaces are two-dimensional except for the eigenspace with eigenvalue 
A = 0, which is one-dimensional. 

If j 7^ and j ^ n/2, ej and e„_j form a basis for the eigenspace corre- 
sponding to eigenvalue Xj. It is not difficult to verify that 

I 1 \ 

cos(7rj'/n) 

cos(27rj/n) 
ycos((n — l)irj/n) J 

I ° \ 

sin(7rj/n) 

sin(27r//n) 
\sin((n- l)nj/n)J 



and 



- (e^ + e„_ j ) 



•■n-jj 
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form a real basis for this cigcnspacc. Let 

ojj = \f— \j = 2sm(nj/n). 

In the case where n is odd, we can write the general solution to our dynamical 
system (2.17) as 

(n-l)/2 

x = e (a + bot) + ^ - (ej + e n ^j ) (aj cos Ujt + bj sin ujjt) 
3=1 

(n-l)/2 



The motion of the system can be described as a superposition of several modes 
of oscillation, the frequencies of oscillation being 



2ir 



k sin(7rj/n) 



m 7r 



Note that the component 



e (a + bot) 



AN 
i 

i 

W 



(a + b t) 



corresponds to a constant speed motion of the carts around the track. If n is 
large, 

sin(7r/n) . ir/n 1 uj\ . Ik 1 

7T TT Tl llT V Tfi Ti ' 

so if we were to set fc/m = n 2 , the lowest nonzero frequency of oscillation would 
approach one as n — > oo. 

Exercise: 

2.5.1. Find the eigenvalues of the matrix 

/0 1 
1 






\1 



0\ 



1 

0/ 
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Figure 2.8: A linear array of carts and springs. 



by expanding the determinant 

-A 1 ••• 
-A 1 ••• 
-A ••• 

••• 1 

1 ••■ -A 

2.6 A linear array of weights and springs* 

Suppose more generally that a system of n — 1 carts containing identical weights 
of mass to, and connected by identical springs of spring constant k, are moving 
along a friction-free track, as shown in Figure 2.7. Just as in the preceding 
section, we can show that the carts will move in accordance with the linear 
system of differential equations 

d 2 x k 
— = A-x = -Fx, 
at z m 

where 





( 


-2 


1 • • 


' \ 
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-2 1 
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p = 







1 -2 •■ 











■• 


"2/ 



We take k = n and to = (1/n), so that k/m = n 2 . 

We can use the following Mathematica program to find the eigenvalues of 
the n x n matrix A, when n = 6: 

n = 6 ; 

m = Table [Max [2-Abs [i-j] ,0] , {i,n-l} ,{j,n-l} ]; 
p = m - 4 IdentityMatrix [n-1] ; 
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a = nA2 p 

eval = Eigenvalues [N [a] ] 

Since 

'=2 i£j = i, 

2-\i-j\i = l i£j = i±l 

<0 otherwise, 

the first line of the program generates the (n — 1) x (n — l)-matrix 
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The next two lines generate the matrices P and A 



-P. Finally, the last 



line finds the eigenvalues of A. If you run this program, you will note that all 
the eigenvalues are negative and that Mathematica provides the eigenvalues in 
increasing order. 

We can also modify the program to find the eigenvalues of the matrix when 
n = 12, n = 26, or n — 60 by simply replacing 6 in the top line with the new 
value of n. We can further modify the program by asking for only the smallest 
eigenvalue lambda [ [n] ] and a plot of the corresponding eigenvector: 

n = 14; 

m = Table [Max [2-Abs [i-j] ,0] , {i,n-l} ,{j,n-l} ]; 

p = m - 4 IdentityMatrix [n-1] ; a = nA2 p; 

eval = Eigenvalues [N [a] ] ; evec = Eigenvectors [N [a] ] ; 

Print [eval [[n-1]]] ; 

ListPlot [evec [ [n-1] ] ] ; 

If we experiment with this program using different values for n, we will see 
that as n gets larger, the smallest eigenvalue seems to approach — n 2 and the 
plot of the smallest eigenvector looks more and more like a sine curve. Thus the 
fundamental frequency of the mechanical system seems to approach w/2ir = 1/2 
and the oscillation in the fundamental mode appears more and more sinusoidal 
in shape. 

When n is large, we can consider this array of springs and weights as a model 
for a string situated along the x-axis and free to stretch to the right and the left 
along the rc-axis. The track on which the cart runs restricts the motion of the 
weights to be only in the x-direction. A more realistic model would allow the 
carts to move in all three directions. This would require new variables yi and 
Zi such that 

Ui(t) = the y-component of displacement of the i-th weight, 
Zi(t) = the z-component of displacement of the z-th weight. 
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Figure 2.9: Shape of the fundamental mode (n=100). 



If we were to set 













(zi\ 












Z2 


X = 




, y = 




, z = 






\x n J 




\x n J 




\Z n J 



then an argument similar to the one given above would show that the vectors 
x, y and z would satisfy the matrix differential equations 

dx k A dy k . dz k . 
— = —Ax, -f- = —Ay, — = —Az. 
dt m dt m dt m 

and each of these could be solved just as before. 

Using a similar approach, we could consider even more complicated systems 
consisting of many masses and springs connected in a two- or three-dimensional 
array. 
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Chapter 3 

Fourier Series 



3.1 Fourier series 

The theory of Fourier series and the Fourier transform is concerned with dividing 
a function into a superposition of sines and cosines, its components of various 
frequencies. It is a crucial tool for understanding waves, including water waves, 
sound waves and light waves. Suppose, for example, that the function f(t) 
represents the amplitude of a light wave arriving from a distant galaxy. The 
light is a superposition of many frequencies which encode information regarding 
the material which makes up the stars of the galaxy, the speed with which the 
galaxy is receding from the earth, its speed of rotation, and so forth. Much of 
our knowledge of the universe is derived from analyzing the spectra of stars and 
galaxies. Just as a prism or a spectrograph is an experimental apparatus for 
dividing light into its components of various frequencies, so Fourier analysis is 
a mathematical technique which enables us to decompose an arbitrary function 
into a superposition of oscillations. 

In the following chapter, we will describe how the theory of Fourier series can 
be used to analyze the flow of heat in a bar and the motion of a vibrating string. 
Indeed, Joseph Fourier's original investigations which led to the theory of Fourier 
series were motivated by an attempt to understand heat flow. 1 Nowadays, the 
notion of dividing a function into its components with respect to an appropriate 
"orthonormal basis of functions" is one of the key ideas of applied mathematics, 
useful not only as a tool for solving partial differential equations, as we will sec 
in the next two chapters, but for many other purposes as well. For example, a 
black and white photograph could be represented by a function f(x, y) of two 
variables, f(x,y) representing the darkness at the point (x,y). The photograph 
can be stored efficiently by determining the components of f(x, y) with respect 
to a well-chosen "wavelet basis." This idea is the key to image compression, 
which can be used to send pictures quickly over the internet. 2 

Courier's research was published in his Theorie analytique de la chaleur in 1822. 

2 See Stephane Mallat, A wavelet tour of signal processing, Academic Press, Boston, 1998. 
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We turn now to the basics of Fourier analysis in its simplest context. A 
function / : R — > R is said to be periodic of period T if it satisfies the relation 

f(t + T) = f(t), for allies. 

Thus f(t) = s'mt is periodic of period 2it. 

Given an arbitrary period T, it is easy to construct examples of functions 
which are periodic of period T — indeed, the function f(t) = sin(^p) is periodic 
of period T because 

. ,2n(t + T). . ,2-Kt „ , . ,27rt. 
sm( \ ' ) = sin(— + 2tt) = sm( — ). 

More generally, if k is any positive integer, the functions 

,2nkt. . ,2irkt. 

cos(-^-) and sm( — ) 

are also periodic functions of period T. 

The main theorem from the theory of Fourier series states that any "well- 
behaved" periodic function of period T can be expressed as a superposition of 
sines and cosines: 

a /27rf /47rf . ,2irt . Ant 

f(t) = — + a lC os( — ) + a 2 cos( — ) + ... + b 1 sm( — ) + b 2 sm( — ) + ... . 

(3.1) 

In this formula, the a^'s and fr^'s are called the Fourier coefficients of /, and 
the infinite series on the right-hand side is called the Fourier series of /. 

Our first goal is to determine how to calculate these Fourier coefficients. For 
simplicity, we will restrict our attention to the case where the period T = 2n, 
so that 

/(*) = ~^ + fl i cos t + a 2 cos 2t + . . . + bi sin t + b 2 sin 2t + . . . . (3.2) 

The formulae for a general period T are only a little more complicated, and are 
based upon exactly the same ideas. 

The coefficient ao is particularly easy to evaluate. We simply integrate both 
sides of (3.2) from — ir to n: 



f(t)dt= / — dt-\- I a\costdt+ I a 2 cos2tdt 

-TT J —TT ^ J— 7T J — TT 

bi sin tdt + / 

-7T J — 7T 



+ 



+ / bi sin tdt + b 2 sin 2tdt + ... . 

J —7T J —IT 

Since the integral of cos kt or sin kt over the interval from — tt to ir vanishes, we 
conclude that 



f{t)dt = 7ra , 
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and we can solve for a , obtaining 



a = - I f(t)dt. (3.3) 



7T 



To find the other Fourier coefficients, we will need some integral formulae. 
We claim that if m and n are positive integers, 

/ 1 7r, for m = n, , . 

cos nt cos ratdt = < (3.4) 
. n |^0, for n, 

. ■ , ( 7r, for m = n, 
sin sin mfai = < (3.5) 
-7T 0, for m ^ n, 



/7T 
-7T 



sin nt cos vntdt = 0. (3.6) 



Let us verify the first of these equations. We will use the trigonometric 
identities 

cos((n + m)t) = cos nt cos mt — sin nt sin mt, 
cos((n — m)t) = cos nt cos mi + sin nt sin mt. 
Adding these together, we obtain 

cos((n + m)t) + cos((n — m)t) = 2 cos nt cos mt, 

or 

cos nt cos mt = -(cos((n + m)t) + cos((n — m)t)). 

Hence 

(cos((n + m)t) + cos((n — m)t))dt, 



/7T 1 Z" 71 

cos nt cos mtdt = - / 
-7T 2 j_. 



and since 

= 0, 



1 

/ cos(kt)dt — — sin(fct) 

J -IT K 



if fc ^ 0, we conclude that 

r 

I cos nt cos mtdt 

J — 7T 



7T, for m = n, 
0, for m 7^ n. 



The reader is asked to verify the other two integral formulae (3.5) and (3.6) 
in the exercises at the end of this section. 
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To find the formula for the Fourier coefficients for k > 0, we multiply 
both sides of (3.2) by sinfct and integrate from — tt to tt: 

/TT P'K ^ f-TT 

f(t) cos ktdt = / — cos ktdt+ / ai cos t cos ktdt +.. . 
-TT J — TT ^ J — IT 



+ . . . I a,k cos kt cos ktdt + . . . 
+ ...+ / bj sin jt cos ktdt + . . . . 

J — TT 

According to our formulae (3.4) and (3.6), there is only one term which survives: 

r 

f(t) cos ktdt = I afc cos kt cos ktdt — na/, ■ 

J — TT 

We can easily solve for a^: 

i r 

a k = - f(t) cos ktdt. (3.7) 

J-TT 

A very similar argument yields the formula 



J-tt 



sin ktdt. 



(3.6 



Example. Let us use formulae (3.3), (3.7), and (3.8) to find the Fourier coeffi- 
cients of the function 



/(*) = 



— TT, for — TT < t < 0, 
tt, for < t < TT, 

0, for t = 0, tt, 



extended to be periodic of period 2tt. Then 



and 



f /(*) 

J —TT 



ao _ 
2 

cos mtdt 

1 -TT 



average value of / = 0, 





1 




tt cos mtdt 


+ - 


/ tt cos mtdt 




TT 


Jo 



1 TT 



[smmt]_ n + - — [sin mt]5 
tt m tt m 



while 



b m = 



f(t) sin mtdt 



-7rsin mtdt 



f 

Jo 



tt sin mtdt 
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Figure 3.1: A graph of the Fourier approximation <p5(t) — 4sint+ (4/3) sin3i + 
(4/5) sin 5i. 



1 7T n 1 — 7T 2 2 

1 cos mt\_- H cosmtL = cosmir 

7r m 7r to mm 

= 1(1 _(_!)"») = /m' if™ is odd, 
to I 0, if to is even. 



Thus we find the Fourier series for /: 



4 4 

/(f) = 4 sin t + - sin 3t + - sin 5t H . 

3 5 

The trigonometric polynomials 

4 

(j>i (t) = 4 sin t, 4>z(t) — Asmt + - sin3t 

o 

4 4 
<fc (*) = 4 sin t + - sin 3i + - sin 5t 
3 5 

are approximations to f(t) which improve as the number of terms increases. 

By the way, this Fourier series yields a curious formula for tt. If we set 
x = ir/2, f(x) — n, and we obtain 

4 4 

tt = 4sin(7r/2) + - sin(37r/2) + - sin(57r/2) H 

3 5 

from which we conclude that 

1111 

- = 4 C 1 -3 + 5-7 + 9--)" 



GG 



-w\/ 



3 



WW 



2 



1 



-4 



-2 



2 



4 



-1 



-2 





Figure 3.2: A graph of the Fourier approximation <fii3 (t). The overshooting near 
the points of discontinuity is known as the "Gibbs phenomenon." 

The Fourier series on the right-hand side of (3.1) is often conveniently ex- 
pressed in the £ notation, 



just as we did for power series in Chapter 1. 

It is an interesting and difficult problem in harmonic analysis to determine 
how "well-behaved" a periodic function f(t) must be in order to ensure that 
it can be expressed as a superposition of sines and cosines. An easily stated 
theorem, sufficient for many applications is: 

Theorem 1. If f is a continuous periodic function of period T, with condin- 
uous derivative f , then f can be written uniquely as a superposition of sines 
and cosines, in accordance with (3.9), where the a^'s and bk's are constants. 
Moreover, the infinite series on the right-hand side of (3.9) converges to f(t) for 
every choice of t. 

However, often one wants to apply the theory of Fourier series to functions 
which are not quite so well-behaved, in fact to functions that are not even be 
continuous, such as the function in our previous example. A weaker sufficient 
condition for / to possess a Fourier series is that it be piece wise smooth. 

The technical definition goes like this: A function f(t) which is periodic of 
period T is said to be piecewise smooth if it is continuous and has a continuous 
derivative /'(i) except at finitely many points of discontinuity within the interval 
[0,T], and at each point t of discontinuity, the right- and left-handed limits of 




(3.9) 
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/ and /' 



— >to+ 



lim (/(*)), 



lim (f(f)), lim (/'(f)), 

: — »to + t — >Iq — 



all exist. The following theorem, proven in more advanced books, 3 ensures 
that a Fourier decomposition can be found for any function which is piecewise 
smooth: 

Theorem 2. If f is any piecewise smooth periodic function of period T, f 
can be expressed as a Fourier series, 



where the <ik 's and bk 's are constants. Here equality means that the infinite 
sum on the right converges to f(t) for each t at which f is continuous. If / is 
discontinuous at to, its Fourier series at to will converge to the average of the 
right and left hand limits of f as t — > to- 

Exercises: 

3.1.1. The function f(t) = cos 2 t can be regarded as either periodic of period 
7r or periodic of period 2tt. Choose one of the two periods and find the Fourier 
series of /(f) . (Hint: This problem is very easy if you use trigonometric identities 
instead of trying to integrate directly.) 

3.1.2. The function /(f) = sin 3 f is periodic of period 2n. Find its Fourier series. 

3.1.3. The function 



can be extended to be periodic of period 2ir. Find the Fourier series of this 
extension. 

3.1.4. The function 



can be extended to be periodic of period 2ir. Find the Fourier series of this 
extension. 

3.1.5. Find the Fourier series of the following function: 



3.1.6. Establish the formulae (3.5) and (3.6), which were given in the text. 

3 See, for example, Ruel Churchill and James Brown, Fourier series and boundary value 
problems, 4th edition, McGraw-Hill, New York, 1987 or Robert Seeley, An introduction to 
Fourier series and integrals, Benjamin, New York, 1966. 



,, . cio v— > 2nkt , . ,2nkt. 

/(*) = y + 2^ a fc c °s(^H + 2^ bkSm (^f^> 
fc=i ' fc=i 




f, for — 7T < f < 7T, 

0, for t = 7r, 



fit) = \t\, 



for f e 
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3.2 Inner products 



There is a convenient way of remembering the formulae for the Fourier coeffi- 
cients that we derived in the preceding section. Let V be the set of piecewise 
smooth functions which are periodic of period 2ir. We say that V is a vec- 
tor space because elements of V can be added and multiplied by scalars, these 
operations satisfying the same rules as those for addition of ordinary vectors 
and multiplication of ordinary vectors by scalars. We define an "inner product" 
between elements of V by means of the formula 



</,<?> 



- r f(t)g(t)dt. 

K J -TV 



Thus for example, if f(t) = s'mt and g(t) = 2 cost, then 

(f,g)=-[ 2 sin t cos tdt = [ sin(2t)dt = - cos(2t)|!^ = 0. 

The remarkable fact is that this inner product has properties quite similar 
to those of the standard dot product on R n : 

• (/) 9) = (flS /)i whenever / and g are elements of V. 

• (f + g,h) = {f,h) + {g,h). 

• ( c /i 9) — c (/) g)> when c is a real constant. 

• (/j /) > 0, with equality holding only if / = (at all points of continuity). 

This suggests that we might use geometric terminology for elements of V 
just as we did for vectors in R n . Thus, for example, we will say that an element 
/ of V is of unit length if (/, /) = 1 and that two elements / and g of V are 
perpendicular if (f,g) = 0. 

In this terminology, the formulae (3.4), (3.5), and (3.6) can be expressed by 
stating that the functions 



l/v2, cost, cos2t, cos3t, ... , sint, sin2t, sin3t, 

are of unit length and perpendicular to each other. Moreover, by the theorem in 
the preceding section, any element of of V can be written as a (possibly infinite) 
superposition of these functions. We will say that the functions 

eo(t) = —i=, ei(t) = cost, e 2 (i) = cos2t, 
v2 

ei(t)=sini, e 2 (i)=sin2t, 

make up an orthonormal basis for V. 

We saw in Section 2.3 that if bi, b2, ... , b„ is an orthonormal basis for R n , 
then 

f€R n f = (f-bi)bi + ••■ + (f -b„)b n . 
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The same formula holds when R n is replaced by V and the dot product is 
replaced by the inner product (•, •): If / is any element of V, we can write 

f(t) = (/(*), e (t))e (t) + (f(t), ei (t)) ei (t) + (/(*), e 2 (t))e 2 (t) + ■■■ 

+ </(*), ei(t))ei(t) + (f(t),e 2 (t))e 2 (t) + ■■■. 

In other words, 

f(t) = (/, -J=) -J= + (/, cos t) cos t + (/, cos 2t) cos 2t + . . . 

+(/, sin t) sin t + (/, sin 2t) sin 2t + . . . , 
= ^ + ai cos t + a 2 cos 2t + . . . + b\ sin t + b 2 sin 2t + . . . , 

where 

ao = (/,!), a 1 = (f,cost), a 2 = (/,cos2i), 

6i = </, sin i), 6 2 = </,sin2t), .... 

Use of the inner product makes the formulae for Fourier coefficients almost 
impossible to forget. 

We say that a function /(i) is odd if /(— t) = —f(t), even if /(— t) = f(t). 
Thus 

sin£, sin2£, sin3i, 
are odd functions, while 

1, cost, cos2i, cos3i, 

are even functions. Let 

^odd = {/ e F : / is odd}, Weven = {/ e V : / is odd}. 

Then 

/,S G Wodd / + .9 e W odd and c/ e W odd , 

for every choice of real number c. Thus we can say that W odd is a linear 
subspace of V. Similarly, Weven is a linear subspace of V. 
It is not difficult to show that 

f^ W odd> 9 & Weven => (f,g) = 0; 

in other words, the linear subspaces W odd and Weven are orthogonal to each 
other. Indeed, under these conditions, fg is odd and hence 

(/,<?> = " I' f(t)g(t)dt=- [ f(t)g(t)dt+- r f(t)g(t)dt 
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= - f f(-t)g(-t)(-dt) + - [ W f(t)g(t)dt 
n J* t Jo 

=-- r -(f(t)g(m-dt)+- r f(t) 9 (t)dt=o. 

T JO K Jo 

The variable of integration has been changed from t to — t in the first integral 
of the second line. 

It follows that if / G W l ^, 

a = (/, 1) = 0, a n = (/, cosnx) = 0, for n > 0. 

Similarly, if / G Weven, 

b n — (/, sinnx) = 0. 

Thus for an even or odd function, half of the Fourier coefficients are auto- 
matically zero. This simple fact can often simplify the calculation of Fourier 
coefficients. 



Exercises: 

3.2.1. Evaluate the inner product 



(f,g) = - r f(t)g(t)dt, 

K J -TV 

in the case where f(t) — cost and g(t) = \ sint|. 
3.2.2. Evaluate the inner product 



(f,9) = [ -f{x)g(x)dx, 
Ji x 



in the case where f(x) = cos(log(a;)) and g(x) = 1, where log denotes the natural 
or base e logarithm. (Hint: Use the substitution x = e u .) 

3.2.3. Determine which of the following functions are even, which are odd, and 
which arc neither even nor odd: 

a. f(t) = t 3 + it. 

b. f(t)=t 2 + \t\. 

c. f(t) = e*. 

d. /(t) = i(e* + e-*). 
c. /(t) = |(e*-e-*). 

f. f(t) = Jo(t), the Bessel function of the first kind. 
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3.3 Fourier sine and cosine series 

Let / : [0, L] — > R be a piecewise smooth function which vanishes at and L. 
We claim that we can express f(t) as the superposition of sine functions, 

f(t) = h sin(7rf/L) + b 2 sin(27rt/L) + ... + &„ sin(n7rf/L) + . . . . (3.10) 

We could prove this using the theory of even and odd functions. Indeed, we can 
extend / to an odd function / : [— L, L] — > R by setting 



/(*) 



f(t), forte[0,L], 
-/(-*), forte[-L,0], 



then to a function / : R — > R, which is periodic of period 2L by requiring that 

f(t + 2L) = f(t), for all te«. 

The extended function lies in the linear subspace WqJ^. It follows from the 
theorem in Section 3.1 that / possesses a Fourier series expansion, and from the 
fact that / is odd that all of the a„'s are zero. On the interval [0, L], / restricts 
to / and the Fourier expansion of / restricts to an expansion of / of the form 
(3.10) which involves only sine functions. We call (3.10) the Fourier sine series 
off. 

A similar argument can be used to express a piecewise smooth function 
/ : [0, L\ — > R into a superposition of cosine functions, 

/(*) = ~^ + a i cos(7rf/L) + a 2 cos(27rf/£) + ... + a n cos{mrt/L) + . . . . (3.11) 

To obtain this expression, we first extend / to an even function / : [-L, L] — > R 
by setting 

'/(*), forte[0,L], 
/(-*), forte[-L,0], 



/(*) 



then to a function / : R — > R which is periodic of period 2L, by requiring that 

f{t + 2L) = f(t), for all t G RR. 

This time the extended function lies in the linear subspace Weven- It follows 
from the theorem in Section 3.1 that / possesses a Fourier series expansion, and 
from the fact that / is even that all of the &„'s are zero. On the interval [0, L], 
f restricts to / and the Fourier expansion of / restricts to an expansion of / of 
the form (3.11) which involves only cosine functions. We call (3.11) the Fourier 
cosine series of /. 

To generate formulae for the coefficients in the Fourier sine and cosine series 
of /, we begin by defining a slightly different inner product space than the one 
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considered in the preceding section. This time, we let V be the set of piecewise 
smooth functions / : [0, L] — > R and let 

V = {f€V: /(0) = = f(L)}, 

a linear subspace of V. We define an inner product (-, •) on V by means of the 
formula 

2 ^ 



(f,g) = z y o f(t)g(t)dt. 

This restricts to an inner product on VJ> 

Let's consider now the Fourier sine series. We have seen that any element 
of Vo can be represented as a superpostion of the sine functions 

sin(7rf/L), sin(27rf/L), sin(nnt/L), .... 

We claim that these sine functions form an orthonormal basis for Vq with respect 
to the inner product we have defined. Recall the trigonometric formulae that 
we used in §3.1: 

cos((n + m)irt/L) — cos(mrt/L) cos{mirt/L) — sin(n7rf / L) s\n(rmrt/L), 

cos((n — m)irt/L) = cos(nnt/L) cos{mirt / L) + sin(n7rf/L) sm(rmrt/L). 
Subtracting the first of these from the second and dividing by two yields 

sm(nnt/L) sin(m7rf/L) = ^(cos((n — m)irt/L) — cos((n + rnjwt/L), 



and hence 

t-L 



sm(nirt/L) sin(mirt/ L)dt = 



L 

(cos((n — rn)7rt/L) — cos((n + rnjirt/ L)dt. 



If n and m are positive integers, the integral on the right vanishes unless n = m, 
in which case the right-hand side becomes 



1 f , L 

2 7o 2 



Hence 



2 f n . , i t \ • / i t \ i I 1, for m = n, 
— / sm(nirt/L)s\n(rmrt/L)dt = { (3.12) 
L J 0, for m ^ n 



Therefore, just as in the previous section, we can evaluate the coefficients of 
the Fourier sine scries of a function / G Vo by simply projecting / onto each 
element of the orthonormal basis. When we do this, we find that 

f(t) = h sin(7rf/L) + b 2 sin(27rf/L) + . . . + sin(mrt/L) + . . . , 
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where 



K = (/,sin(n7rf/L)) 



L 



f(t) sin(mrt/L)dt. 



(3.13) 



We can treat the Fourier cosine series in a similar fashion. In this case, we 
can show that the functions 



V2' 



cos(nt/L), cos(27rf/L), 



cos(nirt/L), 



form an orthonormal basis for V. Thus we can evaluate the coefficients of the 
Fourier cosine series of a function / € V by projecting / onto each element of 
this orthonormal basis. We will leave it to the reader to carry this out in detail, 
and simply remark that when the dust settles, one obtains the following formula 
for the coefficient a n in the Fourier cosine series: 



2 r 

a n =— I f(t)cos(nirt/L)dt. 
L Jo 

Example. First let us use (3.13) to find the Fourier sine series of 
/(*) 



(3.14) 



t, for < t < tt/2, 

7T - t, for TT/2 < t < 7T. 



(3.15) 



In this case, L = n, and according to our formula, 



K = 



I 



it/2 



t sin ntdt 



/ (tt - t) si 

J-kII 



sin ntdt 



We use integration by parts to obtain 

r 7r/2 

t sin ntdt = 



while 



~-t 




tt/2 . 


— cos nt 




+ / 


n 




o Jo 



i-tv/2 



■ cos ntdt 



_ -7rcos(n7r/2) 1 w/2 _ -tt cos(mr/2) sin(n7r/2) 

~ Yn +^^ n ^lo - Yn ' 



[ (w-t) si 

Jtt/2 



sin ntdt - 



-(ir-t) 



cos(nt) 



r i 

rr/2 Jtt/2 « 



cos ntdi 



7rcos(n7r/2) 1 



2n 



- — [smnt]\l /2 = 



7rcos(n7r/2) sin(n7r/2) 



Thus 



b n = 



^ 2 2n 
4sin(n7r/2) 



+ 
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0.5 1 1.5 2 2.5 3 



Figure 3.3: A graph of the Fourier sine approximations <j>i, (f> 3 , cf) 5 and §7. 



and the Fourier sine series of f(t) is 

4 4 4 4 

fit) = — sin t — — sin 3t + —— sin 5t — — — sin 7t + . . . 
J w 7T 9tt 25tt 49tt 

The trigonometric polynomials 

4 4 4 
<j>i (t) = — sin t , 4>3 (t) = — sin t sin 3t, . . . 

7T 7T 97T 

are better and better approximations to the function f(t). 
To find the Fourier cosine series of 



(3.16) 



/(*) = 



t, 



for < t < 7r/2, 



7T — t, for 7r/2 < t < 7T, 



we first note that 



a 



To find the other a„'s, we use (3.14): 



/ t cos ntdt + (n — t) cos nidi 

JO Jir/2 



This time, integration by parts yields 



rir/2 

io ' 



f cos ntdt 



sin nt 



n 



tt/2 



sin ntdt 
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while 



Thus when n > 1, 



7rsin(ri7r/2) cos(n7r/2) 1 
2n n 2 ri 2 



P7T 

/ (7r — t) cos nicft = 



- sin(ni) 



+ / -si 

tt/2 A/2 n 



sin nidt 



-7rsin(n7r/2) 1 

^ ^[oosn*]!^ 

— 7rsin(n7r/2) cos(n7r/2) 1 



2n 



+ 



cos(ri7r). 



: [2cos(n7r/2) - 1 - l(-l) n ], 



and the Fourier sine series of f(t) is 
/(*) 



2 „ 2 2 

cos 2t cos 6t cos Wt 

4 7T 97T 257T 



(3.17) 



Note that we have expanded exactly the same function f(t) on the interval [0, tt] 
as either a superposition of sines in (3.16) or as a superposition of cosines in 
(3.17). 

Exercises: 

3. 3.1. a. Find the Fourier sine series of the following function defined on the 
interval [0, 7r]: 

' At, for < t < it/2, 

4tt - At, for n/2<t< tt. 



/(*) 



b. Find the Fourier cosine series of the same function. 

3. 3. 2. a. Find the Fourier sine series of the following function defined on the 
interval [0,n]: 

f(t)=t(7T-t). 



b. Find the Fourier cosine scries of the same function. 

3.3.3. Find the Fourier sine series of the following function defined on the 
interval [0, 10]: 

't, for < t < 5, 

10 - t, for 5 < t < 10. 



/(*) = 
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3.3.4. Find the Fourier sine series of the following function defined on the 
interval [0,1]: 

/(*) = 5t(l-t). 

3.3.5. (For students with access to Mathematica) Find numerical approximations 
to the first ten coefficients of the Fourier sine series for the function 

f{t)=t + t 2 -2t\ 

defined for t in the interval [0, 1], by running the following Mathematica program 

f[n_] := 2 NIntegrate[(t + tA2 - 2 tA3) Sin[n Pi t] , -ft, 0,1}]; 
b = Table [f [n] , {n,l,10}] 

3.4 Complex version of Fourier series* 

We have seen that if / : R — > R is a well-behaved function which is periodic of 
period 2tt, f can be expanded in a Fourier series 

/(*) — y + a i cos t + a 2 cos(2t) + . . . 

+&i sint + 6 2 sin(2t) + . . . . 

We say that this is the real form of the Fourier series. It is often convenient to 
recast this Fourier series in complex form by means of the Euler formula, which 
states that 

e lB — cos6» + isin6>. 
It follows from this formula that 

e w + e - lB = 2 cos 9, e~ ie + e~ i9 = 2i sin 0, 

or 

e i6 + e -i6 e t0 + e- ie 

cos = , sin 6 = 



2 2i 
Hence the Fourier expansion of / can be rewritten as 

/(*) = y + y ^ + e ^ + f + e ^ + ■ ■ ■ 

+ |(e«- e -«) + |( e 2 «- e - 2 «) + ... , 

or 

/(*) - . . . + c_ 2 e- 24 * + c_ie-** + c + Cl e 4t + c 2 e 2 " + . . . , (3.18) 
where the c^'s are the complex numbers defined by the formulae 
ao ai — ib\ a 2 — i6 2 
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C-l 



ii + ih\ 



C-2 = 



a 2 + i& 2 



If k ^ 0, wc can solve for a k and 6 fc : 

a k = c k + c- k , b k = i{c k - c_fe). (3.19) 
It is not difficult to check the following integral formula via direct integration: 



/" 



e tmt e -int dt 



2ir, for m = n, 
for m ^ n. 



(3.20) 



If we multiply both sides of (3.18) by e jfct , integrate from — tt to 7r and apply 
(3.20), we obtain 



J f(t)e- M dt = 2Trc k , 

J —IT 



which yields the formula for coefficients of the complex form of the Fourier 
series: 



i r 



- M dt. 



(3.21) 



Example. Let us use formula (3.21) to find the complex Fourier coefficients of 
the function 

/(*) = t for -7T < t < 7T, 

extended to be periodic of period 27r. Then 



i r 7 ' 

2^7-. 



te- Lkt dt. 



We apply integration by parts with u = t, dv = e lkt dt, du = dt and w = 
)e" 

1 



(i/k)e- lku 



Ck = 



2tt 



(it/k)e 



-ikt I7T 



(z/fc)e- Jfct * 



1 

27 



7TI 



-ink 



TTl 



ink 



It follows from (3.19) that 



a k = (c k + c_ fe ) = 0, b k = i(c k - c_ fe ) = -2 



(-11 
k 



It is often the case that the complex form of the Fourier series is far simpler to 
calculate than the real form. One can then use (3.19) to find the real form of 
the Fourier series. 

Exercise: 

3.4.1. Prove the integral formula (3.20), presented in the text. 



78 



3.4.2. a. Find the complex Fourier coefficients of the function 

f(t) = t 2 for -7T < t < 7T, 

extended to be periodic of period 2n. 

b. Use (3.19) to find the real form of the Fourier series. 

3.4.3. a. Find the complex Fourier coefficients of the function 

f(t) = t(n - t) for -7T < t < ir, 
extended to be periodic of period 2ir. 
b. Use (3.19) to find the real form of the Fourier series. 

3.4.4. Show that if a function / : R — > R is smooth and periodic of period 2nL, 
we can write the Fourier expansion of / as 

f(t) = ... + c_ 2 e- 2lt / L + c^e- lt ' L + c + Cl e lt / L + c 2 e 2U ' L + ..., 

where 

Cfe 



l — p L f{t)e- M ' L dt. 

* L J-irL 



3.5 Fourier transforms* 

One of the problems with the theory of Fourier series presented in the previous 
sections is that it applies only to periodic functions. There are many times when 
one would like to divide a function which is not periodic into a superposition of 
sines and cosines. The Fourier transform is the tool often used for this purpose. 

The idea behind the Fourier transform is to think of f(t) as vanishing outside 
a very long interval [— 7rL,7rL]. The function can be extended to a periodic 
function f{t) such that f(t + 2irL) = f(t). According to the theory of Fourier 
series in complex form (see Exercise 3.4.4), 

f(t) = ... + c_ 2 e- 2lt ' L + c^e- lt ' L + c + Cl e lt / L + c 2 e 2rf / L + . . . , 

where the c^'s are the complex numbers. 

Definition. If / : R — > R is a piecewise continuous function which vanishes 
outside some finite interval, its Fourier transform is 

/oo 
f{t)e-**dt. (3.22) 
-OO 

The integral in this formula is said to be improper because it is taken from — oo 
to oo; it is best to regard it as a limit, 

/oo pivL 
f(t)e- l ^dt= lim / f(t)e-^dt. 
-oo L ^°°J-nL 
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To explain (3.22), we suppose that / vanishes outside the interval [— irL,irL]. 
We can extend the restriction of / to this interval to a function which is periodic 
of period 2nL. Then 

/oo pirL 
f{t)e- tkt ' L dt= / f{t)e^ kt ' L dt 
-oo J—ttL 

represents 2nL times the Fourier coefficient of this extension of frequency k/L; 
indeed, it follows from (3.21) that we can write 

f(t) = ... + c_ 2 e- 2U / L + c^e- lt ' L + c + Cl e lt ' L + c 2 e 2U ' L + ..., 

for t € [— nL, -jtL], where 



or alternatively, 



/(*) = • • • + ^/(-2/L) e - M / L + -^/(-l /L)e-^ L + 

+ 2^Z /(0) + 2^L /(1/T)e4t/i + 2^L /(2/L)e2 " /L + ' • ' ' 
In the limit as I — > oo, it can be shown that this last sum approaches an 
improper integral, and our formula becomes 

/(*)=2^ j JttY^dt (3.23) 

Equation (3.23) is called the Fourier inversion formula. If we make use of 
Eulcr's formula, we can write the Fourier inversion formula in terms of sines 
and cosines, 

1 f°° - i f°° - 

/(*) = ^ J m) cos ztdt +—J m sm #de, 

a superposition of sines and cosines of various frequencies. 

Equations (3.22) and (3.22) allow one to pass back and forth between a given 
function and its representation as a superposition of oscillations of various fre- 
quencies. Like the Laplace transform, the Fourier transform is often an effective 
tool in finding explicit solutions to differential equations. 

Exercise: 

3.5.1. Find the Fourier transform of the function f(t) defined by 



/(*) 



1, if — 1 < * < 1, 
0, otherwise. 
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Chapter 4 

Partial Differential 
Equations 

4.1 Overview 

A partial differential equation is an equation which contains partial derivatives, 
such as the equation 

du d 2 u 
dt dx 2 ' 

in which u is regarded as a function of x and t. Unlike the theory of ordinary 
differential equations which centers upon one key theorem — the fundamental 
existence and uniqueness theorem — there is no real unified theory of partial dif- 
ferential equations. Instead, each type of partial differential equations exhibits 
its own special features, which usually mirror the physical phenomena which 
the equation was first used to model. 

Many of the foundational theories of physics and engineering are expressed 
by means of systems of partial differential equations. The reader may have 
heard some of these equations mentioned in previous courses in physics. Fluid 
mechanics is often formulated by the Euler equations of motion or the so-called 
Navier-Stokes equations, electricity and magnetism by Maxwell's equations, gen- 
eral relativity by Einstein's field equations. It is therefore important to develop 
techniques that can be used to solve a wide variety of partial differential equa- 
tions. 

In this chapter, we will give two important simple examples of partial dif- 
ferential equations, the heat equation and the wave equation, and we will show 
how to solve them by the techniques of "separation of variables" and Fourier 
analysis. Higher dimensional examples will be given in the following chapter. 
We will see that just as in the case of ordinary differential equations, there is an 
important dichotomy between linear and nonlinear equations. The techniques 
of separation of variables and Fourier analysis are effective only for linear par- 
tial differential equations. Nonlinear partial differential equations are far more 
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difficult to solve, and form a key topic of contemporary mathematical research. 1 
Our first example is the equation governing propagation of heat in a bar of 
length L. We imagine that the bar is located along the x-axis and we let 

u(x, t) = temperature of the bar at the point x at time t. 

Heat in a small segment of a homogeneous bar is proportional to temper- 
ature, the constant of proportionality being determined by the density and 
specific heat of the material making up the bar. More generally, if a(x) denotes 
the specific heat at the point x and p(x) is the density of the bar at x, then the 
heat within the region D Xl}X2 between x\ and x 2 is given by the formula 

Heat within D XlX2 = / p(x)a(x)u(x,t)dx. 

J X\ 

To calculate the rate of change of heat within D Xl:X2 with respect to time, we 
simply differentiate under the integral sign: 



d_ 
dt 



rx 2 1 rx 

/ p(x)a(x)u(x, t)dx = I 

J X\ J J X\ 



■" X2 du 

p{x)a(x) — (x,t)dx. 



Now heat is a form of energy, and conservation of energy implies that if no 
heat is being created or destroyed in D XltX2 , the rate of change of heat within 
D Xl . X2 is simply the rate at which heat enters D XlyX2 . Hence the rate at which 
heat leaves D Xl ^ X2 is given by the exp 

f X2 du 

Rate at which heat leaves D XltX2 = p{x)a(x) — {x,t)dx. (4.1) 

J x x dt 

(More generally, if heat is being created within D Xl _ X2 , say by a chemical reac- 
tion, at the rate p(x)u(x, t) + v(x) per unit volume, then the rate at which heat 



leaves D XuX2 is 



f X2 du f X2 

- / p(x)a(x) — (x,t)dx+ / (p,(x)u(x,t) + v{x))dx.) 

Jx t O'' J x t 

On the other hand, the rate of heat flow F(x, t) is proportional to the partial 
derivative of temperature, 

du 

F{x,t) = -K{x) — {x,t), (4.2) 

where k(x) is the thermal conductivity of the bar at x. Thus we find that the 
rate at which heat leaves the region D XltX2 is also given by the formula 



F(x 2 ,t) -F{x u t) = f 



X2 OF 

— {x,t)dx. (4.3) 



1 Further reading can be found in the many excellent upper-division texts on partial differ- 
ential equations. We especially recommend Mark Pinsky, Partial differential equations and 
boundary-value problems with applications, 2nd edition, McGraw-Hill, 1991. An excellent 
but much more advanced book is Michael Taylor, Partial differential equations: basic theory, 
Springer, New York, 1996. 
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Comparing the two formulae (4.1) and (4.3), we find that 



r 2 dF . , , f X2 . . . . du . 

J ~dx^ X,t ' J P( x ) (T ( x )-g^( x ^) dx - 

This equation is true for all choices of x\ and x 2 , so the integrands on the two 
sides must be equal: 

dF du 
-p{x)a{x)- 



dx dt 



It follows from (4.2) that 



k(x)— — J = — p(a;)<r(a;) - 



dx \ dx J dt 
In this way, we obtain the heat equation 

(<•>£)• <«» 

In the more general case in which heat is being created at the rate p(x)u(x, i) + 
v(x) per unit length, one could show that heat flow is modeled by the equation 

p(x)a(x)^ = ^ {^(^^j + ^ x ) u + v ( x )- ( 4 - 5 ) 

In the special case where the bar is homogeneous , i.e. its properties are the 
same at every point, p(x), a(x) and k(x) are constants, say a and k respectively, 
and (4.4) becomes 

du k d 2 u 
dt pa dx 2 ' 

This is our simplest example of a linear partial differential equation. Although 
its most basic application concerns diffusion of heat, it arises in many other 
contexts as well. For example, a slight modification of the heat equation was 
used by Black and Scholes to price derivatives in financial markets. 2 

Exercises: 

4.1.1. Study of heat flow often leads to "boundary- value problems" for ordi- 
nary differential equations. Indeed, in the "steady-state" case, in which u is 
independent of time, equation (4.5) becomes 

k(x)^-(x) J + (j,(x)u(x) + v(x) = 0, 



dx V dx 



2 A description of the Black-Scholcs technique for pricing puts and calls is given in Paul 
Wilmott, Sam Howison and Jeff Dewynne, The mathematics of financial derivatives, Cam- 
bridge Univ. Press, 1995. 
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a linear ordinary differential equation with variable coefficients. Suppose now 
that the temperature is specified at the two endpoints of the bar, say 

u(0) = a, u(L) = 13. 

Our physical intuition suggests that the steady-state heat equation should have 
a unique solution with these boundary conditions. 

a. Solve the following special case of this boundary- value problem: Find u{x), 
defined for < x < 1 such that 

^ = 0, w(0) = 70, = 50. 

b. Solve the following special case of this boundary- value problem: Find u(x), 
defined for < x < 1 such that 

^-m = 0, w(0) = 70, u(l) = 50. 

c. Solve the following special case of this boundary- value problem: Find u(x), 
defined for < x < 1 such that 

d 2 u 

+ x(l - x) = 0, m(0)=0, u(1) = 0. 

d. (For students with access to Mathematica) Use Mathcmatica to graph the 
solution to the following boundary- value problem: Find u{x), defined for < 
x < 1 such that 

d 2 u 

^2 + (1 + x 2 )u = 0, m(0) = 50, u(1) = 100. 

You can do this by running the Mathcmatica program: 

a = 0; b = 1; alpha = 50; beta = 100; 
sol = NDSolve[ {u" [x] + (1 + xA2) u[x] == 0, 
u[a] == alpha, u[b] == beta.}, u, {x,a,b}]; 
Plot [ Evaluate [ u[x] /. sol], {x,a,b}] 

4. 1.2. a. Show that the function 

u (x,t) = -^e— 2 / 4 ' 

V47Tt 

is a solution to the heat equation (4.6) for t > in the case where n/(pa) = 1. 
b. Use the chain rule with intermediate variables x — x — a, i = t to show that 

u a (x,t) = ^L=e-^ 2 ^ 
V47rf 
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is also a solution to the heat equation. 



c. Show that 

/oo 
uq(x, t)dx = 1. 
-oo 

Hint: Let / denote the integral on the left hand side and note that 
I 2 = — — I / e-^+y^dxdy. 

4?rf J-oo J-oo 

Then transform this last integral to polar coordinates. 

d. Use Mathematica to sketch Uo(x, t) for various values of t. What can you 
say about the behaviour of the function u (x,t) as t — > 0? 

e. By differentiating under the integral sign show that if h : R — ► R is any 
smooth function which vanishes outside a finite interval [—L,L], then 

/oo 
u a (x,t)h(a)da (4.7) 
-oo 

is a solution to the heat equation. 

REMARK: In more advanced courses it is shown that (4.7) approaches h{x) as 
t — > 0. In fact, (4.7) gives a formula (for values of t which are greater than zero) 
for the unique solution to the heat equation on the infinite line which satisfies 
the initial condition u(x, 0) = h(x). In the next section, we will see how to solve 
the initial value problem for a rod of finite length. 



4.2 The initial value problem for the heat equa- 
tion 

We will now describe how to use the Fourier sine series to find the solution to 
an initial value problem for the heat equation in a rod of length L which is 
insulated along the sides, whose ends are kept at zero temperature. We expect 
that there should exist a unique function u(x,t), defined for < x < L and 
t > such that 

1. u(x,t) satisfies the heat equation 

du 9 d 2 u , , „, 

dt= C W> (4 - 8) 



where c is a constant. 

2. u(x,t) satisfies the boundary condition u(0,t) = u(L,t) = 0, in other 
words, the temperature is zero at the endpoints. (This is sometimes called 
the Dirichlet boundary condition.) 
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3. u(x,t) satisfies the initial condition u(x,0) = h(x), where h(x) is a given 
function, the initial temperature of the rod. 



In more advanced courses, it is proven that this initial value problem does in 
fact have a unique solution. We will shortly see how to find that solution. 

Note that the heat equation itself and the boundary condition are homoge- 
neous and linear — this means that if u\ and u 2 satisfy these conditions, so does 
CiUi + c 2 u 2 , for any choice of constants C\ and c 2 . Thus homogeneous linear 
conditions satisfy the principal of superposition. 

Our method makes use of the dichotomy into homogeneous and nonhomo- 
geneous conditions: 

Step I. We find all of the solutions to the homogeneous linear conditions of the 
special form 



By the superposition principle, an arbitrary linear superposition of these solu- 
tions will still be a solution. 

Step II. We find the particular solution which satisfies the nonhomogcneous 
condition by Fourier analysis. 

Let us first carry out Step I. We substitute u(x,t) = f(x)g(t) into the heat 
equation (4.8) and obtain 



Now we separate variables, putting all the functions involving t on the left, all 
the functions involving x on the right: 



The left-hand side of this equation does not depend on x, while the right-hand 
side does not depend on t. Hence neither side can depend upon cither x or t. 
In other words, the two sides must equal a constant, which we denote by A and 
call the separating constant. Our equation now becomes 



u(x,t) = f(x)g(t). 



f(x)g'(t) = <?f"{x)g{t). 



g'(t) _/"(*) 



= A 



c 2 g(t) fix) 



which separates into two ordinary differential equations, 



g'(t) 



A, or g'{t) = Xc 2 g(t) 



(4.9) 



<*g{t) 



and 



f"(x) 



= A, or f"(x) = Xf(x). 



(4.10) 



fix) 
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The homogeneous boundary condition u(0,t) = u(L,t) = becomes 
f(0)g(t) = f(L)g(t) = 0. 
If g(t) is not identically zero, 

/(0) = f(L) = 0. 

(If g(t) is identically zero, then so is u(x, t), and we obtain only the trivial 
solution u = 0.) 

Thus to find the nontrivial solutions to the homogeneous linear part of the 
problem requires us to find the nontrivial solutions to a boundary value problem 
for an ordinary differential equation: 

Hs) = ^j(/(*)) = */(*)> /(°)=° = /( i )- ( 4 - n ) 
We will call (4.11) the eigenvalue problem for the differential operator 




acting on the space Vo of well-behaved functions / : [0, L] — ► R which vanish at 
the endpoints and L. 

We need to consider three cases. (As it turns out, only one of these will 
actually yield nontrivial solutions to our eigenvalue problem.) 

Case 1: A = 0. In this case, the eigenvalue problem (4.11) becomes 

f"(x) = Q, /(0) = = /(L). 

The general solution to the differential equation is f(x) — ax + b, and the only 
particular solution which satisfies the boundary condition is / = 0, the trivial 
solution. 

Case 2: A > 0. In this case, the differential equation 

f"(x) = Xf(x), or f"(x) — Xf(x) = 
has the general solution 

f(x) = cie^ 1 + c 2 e~ VXx . 
It is convenient for us to change basis in the linear space of solutions, using 

cosh(\/A:E) = ^(e^ x + e-^ x ), sin^A/Aa;) = ^(e^ x - e"^) 

instead of the exponentials. Then we can write 

f(x) = acosh(VXx) + bsmh(V\x), 
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with new constants of integration a and 6. We impose the boundary conditions: 
first 

/(0)=0 => a = => f(x) = bsmh(y/\x), 

and then 

f(L) = => 6sinh(\/AL) = => 6 = => f(x) = 0, 
so we obtain no nontrivial solutions in this case. 

Case 3: A < 0. In this case, we set u = \/—\, and rewrite the eigenvalue 
problem as 

f"(x) + Lo 2 f(x) = 0, /(0)=0 = /(L). 

We recognize here our old friend, the differential equation of simple harmonic 
motion. We remember that the differential equation has the general solution 

f(x) — acos(ujx) + 6sin(u;a;). 

Once again 

/(0) = => a = => f(x) = 6sin(wx). 

Now, however, 

f(L) = => 6sin(wL) => b = or sin(wi) = 0, 

and hence either 6 = and we obtain only the trivial solution or sin(wL) = 0. 
The latter possibility will occur if utL = nir, or w = (nw/L), where n is an 
integer. In this case, we obtain 

f(x) = bsm(nirx/L). 

Therefore, we conclude that the only nontrivial solutions to (4.11) are constant 
multiples of 

f[x) = sm(mvx/L), with A = — (mr/L) 2 , n — 1,2,3,.... 

For each of these solutions, we need to find a corresponding g(t) solving equation 
(4.9), 

g'{t) = \c 2 g{t), 

where A = —(nir/L) 2 . This is just the equation of exponential decay, and has 
the general solution 

g(t) = be-i^/L) 2 ^ 

where a is a constant of integration. Thus we find that the nontrivial prod- 
uct solutions to the heat equation together with the homogeneous boundary 
condition u(0,t) — — u(L,t) arc constant multiples of 

u n {x,t) = sm{nnx/L)e- (ncn/L)2t . 



<S<S 



It follows from the principal of superposition that 

u(x, t) = h s\n{Kx/L)e-^l^ 2t + b 2 sm^Ttx/L^-^/^ 1 + ... (4.12) 

is a solution to the heat equation together with its homogeneous boundary 
conditions, for arbitrary choice of the constants b\, b 2 , ... . 

Step II consists of determining the constants b n in (4.12) so that the initial 
condition u(x, 0) = h(x) is satisfied. Setting t = in (4.12) yields 

h(x) = u(x, 0) = 61 sin(7rx/£) + b 2 sin(27rx/i) + . . . . 

It follows from the theory of the Fourier sine series that h can indeed be rep- 
resented as a superposition of sine functions, and we can determine the 6 n 's as 
the coefficients in the Fourier sine series of h. Using the techniques described in 
Section 3.3, we find that 

2 f L 

b n = — h(x)sm(mrx/L)dx. 
L Jo 



Example 1. Suppose that we want to find the function u(x,i), defined for 
< x < 7r and t > 0, which satisfies the initial- value problem: 

— = — u(0, t) = u(n, t) = 0, u(x,0) = h(x) = 4sinx+2sin2.T+7sin3a;. 

In this case, the nonvanishing coefficients for the Fourier sine series of h are 

61 - 4, b 2 = 2, 63 = 7, 

so the solution must be 

u(x, t) = 4sin:re~ t + 2sin2xe~ 4 * + 7sin3a;e~ 9t . 



Example 2. Suppose that we want to find the function u(x,t), defined for 
< x < 7r and t > 0, which satisfies the initial- value problem: 

~di = dx~^' u(0,t) = u(n,t) = 0, u(x, 0) = h(x), 



where 

x, for < x < 7r/2, 

7r — x, for 7r/2 < X < TT. 

We saw in Section 3.3 that the Fourier sine series of h is 



h(x) 



4 4 4 4 

h(x) — — sin x sin 3x H sin 5x sin 7x + 

TT 9tt 25tt 49tt 



and hence 



4 4 4 

u(x, t) — — sin xe~ l — — sin 3xe~ 9t + —— sin 5xe _25t 

7T 9n 257T 
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Exercises: 

4.2.1. Find the function u(x, t), defined for < x < it and t > 0, which satisfies 
the following conditions: 

du d 2 u . 
~Qt = ~Q^2-' u(0,t) = u(w,t) = 0, u(x,0) = sm2x. 

You may assume that the nontrivial solutions to the eigenvalue problem 

/"(x) = A/Ob), /(0)=0 = /(tt) 

are 

A = — n 2 , /(x) = 6„ sin nx, for n = 1, 2, 3, . . . , 
where the &„'s are constants. 

4.2.2. Find the function u(x, t), defined for < x < it and t > 0, which satisfies 
the following conditions: 

du d 2 u 

— = 7—7, u(0, t) — u(tt, t) — 0, mix, 0) = sin.x + 3sin2a; — 5sin3x. 

4.2.3. Find the function u(x, t), defined for < x < n and t > 0, which satisfies 
the following conditions: 

~di ^ dx 2 ' u{0,t) = u{iT,t) = 0, u(x,0) = x(tt - x). 
(In this problem you need to find the Fourier sine series of = x(7r — x).) 

4.2.4. Find the function u(x, t), defined for < x < n and t > 0, which satisfies 
the following conditions: 

1 du d 2 u 

2t + 1 ~di = dx 2 ' u ( '*) = "(tT)*) = °> m(x, 0) = sin x + 3 sin 2x. 

4.2.5. Find the function u(x, t), defined for < x < n and t > 0, which satisfies 
the following conditions: 

du 2 u 

(t + 1)"^ = g~ 2 ' u(0, t) = u(7r, t) = 0, m(x,0) = sinx + 3sin2x. 

4. 2. 6. a. Find the function w(x), defined for < x < ir, such that 

_ = 0, tu(0) = 10, w(?r) = 50. 

b. Find the general solution to the following boundary value problem for the 
heat equation: Find the functions u(x,t), defined for < x < ir and t > 0, such 
that 

Ou 3 2 u 

- d ~t=dx 2 ~' u(0,t) = 10, u(7r,t)=50. (4.13) 



90 



(Hint: Let v = u — w, where w is the solution to part a. Determine what 
conditions v must satisfy.) 

c. Find the particular solution to (4.13) which in addition satisfies the initial 
condition 

40 

u(x, 0) = 10 H x + 2 sin x — 5 sin 2x. 

TT 

4. 2. 7. a. Find the eigenvalues and corresponding eigenfunctions for the differen- 
tial operator 

which acts on the space V of well-behaved functions / : [0,7r] — > R which vanish 
at the endpoints and ir by 

Hfl-g + S/. 

b. Find the function u(x,t), defined for < x < n and t > 0, which satisfies 
the following conditions: 

du d 2 u 

— = — -tt + 3u, u(0, t) = u(tt, t) = 0, u(x, 0) = sin x + 3sin2x. 
at ox 1 

4.2.8. The method described in this section can also be used to solve an initial 
value problem for the heat equation in which the Dirichlet boundary condition 
u(0,t) = u(L, t) = is replaced by the Neumann boundary condition 

g(0, t ) = g(L, t )=0. (4.14) 

(Physically, this corresponds to insulated endpoints, from which no heat can 
enter or escape.) In this case separation of variables leads to a slightly different 
eigenvalue problem, which consists of finding the nontrivial solutions to 

f " {x) = & {f{x)) = A/(2:) ' //(0) = ° = 

a. Solve this eigenvalue problem. (Hint: The solution should involve cosines 
instead of sines.) 

b. Find the general solution to the heat equation 

du d 2 u 
dt dx 2 

subject to the Neumann boundary condition (4.14). 
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c. Find the function u(x, t), defined for < x < n and t > 0, such that: 

~di = d~x 2 ~ 1 dx^ 0,t ^ = dx^'^ = °' u ( x >°) ^Z^x + 'tc^^- 
i^.d. We can also treat a mixture of Dirichlct and Neumann conditions, say 

3u 

u(0,t) = —(L,t) = 0. (4.15) 

In this case separation of variables leads to the eigenvalue problem which consists 
of finding the nontrivial solutions to 

f " {x) = &^ f{x)) = A/(2:) ' /(0) = ° = f{L) - 

a. Solve this eigenvalue problem. 

b. Find the general solution to the heat equation 

du d 2 u 
dt dx 2 

subject to the mixed boundary condition (4.15). 

c. Find the function u(x, t), defined for < x < n and t > 0, such that: 

— =q^2, u{0,t) = —{w,t)=0, u(x,0) =4sin(a;/2) + 12sin(3a;/2). 

4.3 Numerical solutions to the heat equation 

There is another method which is sometimes used to treat the initial value 
problem described in the preceding section, a numerical method based upon 
"finite differences." Although it yields only approximate solutions, it can be 
applied in some cases with variable coefficients when it would be impossible to 
apply Fourier analysis in terms of sines and cosines. However, for simplicity, 
we will describe only the case where p and k are constant, and in fact we will 
assume that c 2 = L = 1. 

Thus we seek the function u(x,t), defined for < x < 1 and t > 0, which 
solves the heat equation 

du d 2 u 
dt dx 2 

subject to the boundary conditions u(0,t) = u(l,t) = and the initial con- 
dition u(x, 0) = h(x), where h(x) is a given function, representing the initial 
temperature. 

For any fixed choice of t the function u(x, to) is an element of Vo, the space of 
piece wise smooth functions defined for < x < 1 which vanish at the endpoints. 
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Our idea is to replace the "infinite-dimensional" space Vo by a finite-dimensional 
Euclidean space J?" -1 and reduce the partial differential equation to a system 
of ordinary differential equations. This corresponds to utilizing a discrete model 
for heat flow rather than a continuous one. 
For < i < n, let Xj = i/n and 

Ui(t) = u(xi,t) = the temperature at Xi at time t. 

Since Uo(t) = = u n (t) by the boundary conditions, the temperature at time t 
is specified by 

/ «i(i) \ 
u(t) = U2(t) 

\U n -i(t)J 

a vector-valued function of one variable. The initial condition becomes 



u(0) 



where 



h = 



f h( Xl ) \ 
h(x 2 ) 



\h(x n -i)J 
derivativ 

du ( ' x t + Xi+i j\ ^ u l+ i(t) - Ui(t) _ [u i+ i(t) - Ui(t)\ 
dx 



We can approximate the first-order partial derivative by a difference quo- 
tient: 



n[u i+ i(t) - Ui(t)]. 



.* = 



Xi+i — Xi l/n 
Similarly we can approximate the second-order partial derivative: 

/ Xj+x i+1 ,\ du ( Xj-l+Xj ,\ 
\ 2 i l ) dx \ 2 ' L J 



dx 2 



(Xi,t) 



du 
dx 



= n 



du f Xi + x 



dx 



l/n 

du ( Xi-i + Xi 
~ da 



= n 2 [u,_!(t) - 2ui(t) + Ui+i(t)]. 
Thus the partial differential equation 

du d 2 u 
dt dx 2 

can be approximated by a system of ordinary differential equations 

dui 
~dt 

This is a first order linear system which can be presented in vector form as 



du 

~dt 



n 2 (ui-i - 2ui + 
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1 • • 
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the last matrix having n — 1 rows and n — 1 columns. Finally, we can rewrite 
this as 

^ = Au, where A = n 2 P, (4.16) 

a system exactly like the ones studied in Section 2.5. In the limit as n — > oo 
one can use the Mathematica program of §2.6 to check that the eigenvalues of A 
approach the eigenvalues of d 2 /dx 2 as determined in the preceding section, and 
the eigenvectors approximate more and more closely the standard orthonormal 
basis of sine functions. 

One could continue constructing a numerical method for solution of our ini- 
tial value problem by means of another discretization, this time in the time 
direction. We could do this via the familiar Cauchy-Euler method for finding 
numerical solutions to the linear system (4.16). This method for finding ap- 
proximate solutions to the heat equation is often called the method of finite 
differences. With sufficient effort, one could construct a computer program, 
using Mathematica or some other software package, to implement it. 

More advanced courses on numerical analysis often treat the finite difference 
method in detail. 3 For us, however, the main point of the method of finite 
differences is that it provides considerable insight into the theory behind the 
heat equation. It shows that the heat equation can be thought of as arising 
from a system of ordinary differential equations when the number of dependent 
variables goes to infinity. It is sometimes the case that either a partial differen- 
tial equation or a system of ordinary differential equations with a large or even 
infinite number of unknowns can give an effective model for the same physical 
phenomenon. This partially explains, for example, why quantum mechanics 
possesses two superficially different formulations, via Schrodinger's partial dif- 
ferential equation or via "infinite matrices" in Hciscnberg's "matrix mechanics." 



4.4 The vibrating string 

Our next goal is to derive the equation which governs the motion of a vibrating 
string. We consider a string of length L stretched out along the x-axis, one end 
of the string being at x = and the other being at x — L. We assume that the 
string is free to move only in the vertical direction. Let 

u(x, t) — vertical displacement of the string at the point x at time t. 

We will derive a partial differential equation for u(x,t). Note that since the 
ends of the string are fixed, we must have u(0, t) = = u(L, t) for all t. 

It will be convenient to use the "configuration space" Vo described in Sec- 
tion 3.3. An element u(x) £ Vo represents a configuration of the string at some 

3 For further discussion of this method one can refer to numerical analysis books, such as 
Burden and Faircs, Numerical analysis, Seventh edition, Brooks Cole Publishing Company, 
2000. 
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instant of time. We will assume that the potential energy in the string when it 
is in the configuration u{x) is 



(4.17) 



where T is a constant, called the tension of the string. 

Indeed, we could imagine that we have devised an experiment that measures 
the potential energy in the string in various configurations, and has determined 
that (4.17) does indeed represent the total potential energy in the string. On 
the other hand, this expression for potential energy is quite plausible for the 
following reason: We could imagine first that the amount of energy in the string 
should be proportional to the amount of stretching of the string, or in other 
words, proportional to the length of the string. From vector calculus, we know 
that the length of the curve u — u(x) is given by the formula 

Length = / \J\ + (du/dx) 2 dx. 
Jo 

But when du/dx is small, 



1 / du 



2 \dx 



^ ) + a small error, 
dx J 



and hence 



\J\ + (du/dx) 2 is closely approximated by 1 + ^(du/dx) 2 . 

Thus to a first order of approximation, the amount of energy in the string should 
be proportional to 



Jo 



_ 1 f du 
1+ 2U 



dx 



-I 



L 1 fduV , 

I — dx + constant. 



2 \dx J 



Letting T denote the constant of proportionality yields 



energy in string = 



L T fdm^ 2 
dx 



dx + constant. 



Potential energy is only defined up to addition of a constant, so we can drop 
the constant term to obtain (4.17). 

The force acting on a portion of the string when it is in the configuration 
u(x) is determined by an element F(x) of Vo- We imagine that the force acting 
on the portion of the string from x to x + dx is F(x)dx. When the force 
pushes the string through an infinitesimal displacement £(x) e Vq, the total 
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work performed by F(x) is then the "sum" of the forces acting on the tiny 
pieces of the string, in other words, the work is the "inner product" of F and £, 

L 



(F(xU(x))= f F(x)£(x)dx. 
Jo 



(Note that the inner product we use here differs from the one used in Section 3.3 
by a constant factor.) 

On the other hand this work is the amount of potential energy lost when the 
string undergoes the displacement: 



2 \dx J J 2 \ dx 

L dudf f L TfdC 2 



T J dxdx dx + J 2 \dx) 

We are imagining that the displacement £ is infinitcsimally small, so terms 
containing the square of £ or the square of a derivative of £ can be ignored, and 
hence 

rL du d£ 



= - T J Q 3xtx dX - 



Integration by parts yields 



L d 2 u „ , N , „ / du A , , N „ / d 



(F(xU(x)) = T J q - i{ x )d x -T^-J) (L) -T^-J) (0). 

Since £(0) = £(L) = 0, we conclude that 

f L f L d 2 u 

^ F(x)^x)dx=(F(xU(x))=Tj Q^i{x)dx. 

Since this formula holds for all infinitesimal displacements we must have 

for the force density per unit length. 

Now we apply Newton's second law, force = mass x acceleration, to the 
function u(x,t). The force acting on a tiny piece of the string of length dx is 
F(x)dx, while the mass of this piece of string is just pdx, where p is the density 
of the string. Thus Newton's law becomes 

J 2 «, , d 2 u 
T _ dx = pdx _. 

If we divide by pdx, we obtain the wave equation, 

d 2 u _ TcPu <Pu _ 2 8 2 u 

dt 2 p dx 2 ' dt 2 dx 2 ' 
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where c 2 = T / p. 

Just as in the preceding section, we could approximate this partial differential 
equation by a system of ordinary differential equations. Assume that the string 
has length L = 1 and set Xi = i/n and 

Ui(t) = u(xi,t) = the displacement of the string at Xi at time t. 

Then the function u(x,t) can be approximated by the vector-valued function 



u(t) 



/ Ui(t) \ 

u 2 {t) 
\u n -i(t)J 



of one variable, just as before. The wave equation is then approximated by the 
system of ordinary differential equations 

d2u 2 2 D 

where P is the (n — 1) x (n — 1) matrix described in the preceding section. Thus 
the differential operator 

L = — g is approximated by the symmetric matrix n 2 P 1 

and we expect solutions to the wave equation to behave like solutions to a 
mechanical system of weights and springs with a large number of degrees of 
freedom. 

Exercises: 

4. 4.1. a. Show that if / : R — > R is any well-behaved function of one variable, 

u(x, t) = f(x + ct) 
is a solution to the partial differential equation 

du du 
dt dx 

(Hint: Use the "chain rule.") 

b. Show that if g : R — > R is any well-behaved function of one variable, 

u(x, t) = g(x — ct) 
is a solution to the partial differential equation 

du du 
dt dx 
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c. Show that for any choice of well-behaved functions / and g, the function 

u(x, t) = f(x + ct) + g(x — ct) 
is a solution to the differential equation 



' d 


d ' 


( du 


du\ 


dt~ 


dx 


\di 


~ C Wx) 



= 0. 



d 2 u 2 d 2 u 
dt 2 dx 2 

Remark: This gives a very explicit general solution to the equation for the 
vibrations of an infinitely long string. 

d. Show that 

f{x + ct) + f{x-ct) 
u(x,t) = 

is a solution to the initial value problem 

d 2 u 9 d 2 u „ . . ,. . du . „. 

W~ c dx 2 ~ = ' u ( x >°) = f( x )> ^^ ) = - 

4.4.2. Show that if the tension and density of a string are given by variable 
functions T(x) and p(x) respectively, then the motion of the string is governed 
by the equation 

d 2 u 1 d f r ,du s 



Of 2 p{x)dx\ T ^'dx 

4.5 The initial value problem for the vibrating 
string 

The Fourier sine series can also be used to find the solution to an initial value 
problem for the vibrating string with fixed endpoints at x — and x = L. 
We formulate this problem as follows: we seek a function u(x,t), defined for 
< x < L and t > such that 

1. u(x,t) satisfies the wave equation 

^ = c 2 ^ (4 18) 

dt 2 dx 2 ' { ' 

where c is a constant. 

2. u(x,t) satisfies the boundary condition u(0,t) = u(L,t) = 0, i.e. the 
displacement of the string is zero at the endpoints. 

3. u(x, t) satisfies the initial conditions 

u(x, 0) = h\(x) and —(x,0) = h 2 (x), 

where h\(x) and h,2(x) are given functions, the initial position and velocity 
of the string. 



<)<S 



Note that the wave equation itself and the boundary condition are homogeneous 
and linear, and therefore satisfy the principal of superposition. 
Once again, we find the solution to our problem in two steps: 

Step I. We find all of the solutions to the homogeneous linear conditions of the 
special form 

u(x,t) = f(x)g(t). 

Step II. We find the superposition of these solution which satisfies the nonho- 
mogcneous initial conditions by means of Fourier analysis. 

To carry out Step I, we substitute u(x,t) = f(x)g(t) into the wave equation 
(4.18) and obtain 

f(x)g"(t) = c 2 f"(x)g(t). 

We separate variables, putting all the functions involving t on the left, all the 
functions involving x on the right: 

g"(t) 2 f"(x) 

g(t) m ■ 

Once again, the left-hand side of this equation does not depend on x, while the 
right-hand side does not depend on t, so neither side can depend upon either x 
or t. Therefore the two sides must equal a constant A, and our equation becomes 

g"(t) /"(*) = Aj 



c 2 g(t) f(x) 
which separates into two ordinary differential equations, 

g"(t) 



c 2 g(t) 

and 



A, or g"{t) = Xc 2 g(t), (4.19) 



= A, or f"(x) = \f(x). (4.20) 



Just as in the case of the heat equation, the homogeneous boundary condition 
u(0,t) = u(L, t) = becomes 

/(0) fl (t) = f(L)g(t) = 0, 

and assuming that g(t) is not identically zero, we obtain 

/(0) = f(L) = 0. 

Thus once again we need to find the nontrivial solutions to the boundary value 
problem, 
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and just as before, we find that the the only nontrivial solutions are constant 
multiples of 

f(x) = sm(nirx/L), with A = — (nir/L) 2 , n = l,2,3, .... 

For each of these solutions, we need to find a corresponding g(t) solving 
equation (4.19), 

g "(t) = -(nir/L) 2 c 2 g(t), or g"(t) + (nn/L) 2 c 2 g(t) = 0. 

This is just the equation of simple harmonic motion, and has the general solution 

g(t) = a cos (ncnt / L) + b sm(ncirt / L) , 

where a and b are constants of integration. Thus we find that the nontrivial 
product solutions to the wave equation together with the homogeneous bound- 
ary condition u(Q, t) = = u(L, t) are constant multiples of 

u n (x,t) — [a n cos(ncnt/L) + b n sin(ncTrt/L)]s\n(mrx/L). 

The general solution to the wave equation together with this boundary condition 
is an arbitrary superposition of these product solutions: 

u(x, t) = [a\ cos(c7rf / L) + bi sin(c7rf/L)] sin(7rx/i) (4-21) 

+ [a 2 cos(2c7ri/L) + b 2 sm(2cirt/L)} sm{2irx/L) + ... . (4.22) 

The vibration of the string is a superposition of a fundamental mode which has 
frequency 

C7T 1 C \JT / p 

~L2^ ~ 2L ~ 2L ' 
and higher modes which have frequencies which are exact integer multiples of 
this frequency. 

Step II consists of determining the constants a n and b n in (4.22) so that the 
initial conditions 

du 

u(x, 0) = h\(x) and — (x, 0) = h,2{x) 

are satisfied. Setting t = in (4.22) yields 

hi(x) = u(x, 0) = ai sin(7ra;/L) + a 2 sin(27rx/L) + . . . , 

so we see that the a„'s are the coefficients in the Fourier sine series of hi. 
If we differentiate equation(4.22) with respect to t, we find that 

Oil . . . C7T . . , C7T . / j- M • / /t\ 

— (x, t) = [ — —ai sm{cirt/L) + —b\ cos(c7rf/_L)J sin(7rx/L) 
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2c7T C7T 

[a 2 s'm(2cirt/L) + — b 2 sm(2cirt/L)] cos(2irx / L) + . . . 



L " L 

and setting t — yields 

h 2 (x) = -wr(x, 0) = —hi sin(7rx/i) H — — b 2 sin(2-7rx/i) + . . . . 
We conclude that 

TIC7T 

—j—b n = the n-th coefficient in the Fourier sine series of h 2 (x). 

Example. Suppose that we want to find the function u(x,t), defined for < 
x < 7r and t > 0, which satisfies the initial- value problem: 

^2=^2' «(0,t)=«(7T,t)=0 ) 

0) = 5 sin x + 12 sin 2a; + 6 sin 3x, — (x, 0) = 0. 
In this case, the first three coefficients for the Fourier sine series of h are 

ai = 5, a 2 = 12, a 3 = 6, 
and all the others are zero, so the solution must be 

u(x, t) = 5 sin x cos t + 12 sin 2x cos 2t + 6 sin 3x cos 3i. 

Exercises: 

4.5.1 What happens to the frequency of the fundamental mode of oscillation of 
a vibrating string when the length of the string is doubled? When the tension 
on the string is doubled? When the density of the string is doubled? 

4.5.2. Find the function u(x, t), defined for < x < n and t > 0, which satisfies 
the following conditions: 

d 2 u d 2 u . , , , 
^2=^2' u(0,t)=u(7r,t)=0, 

u(x, 0) = sin 2x, —(x,0)=0. 
You may assume that the nontrivial solutions to the eigenvalue problem 
f"(x) = Xf(x), /(0)=0 = /(tt) 

are 

A = — n 2 , f(x) = b n sinnx, for n = 1, 2, 3, . . . , 
where the &„'s are constants. 
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4.5.3. Find the function u(x, t), denned for < x < n and t > 0, which satisfies 
the following conditions: 

d 2 u d 2 u , n . . . 
^2=^2' «(0,t)=«(7r,t)=0, 

u(x, 0) = sin a: + 3 sin 2a; — 5 sin 3a;, — (i, 0) = 0. 

at 

4.5.4. Find the function u(x, t), defined for < x < n and t > 0, which satisfies 
the following conditions: 

d 2 u d 2 u in . . . n 
-Q{2=Q^2i «(0,t)=«(7r,t)=0, 

u(x,0) = a;(7r - x), —(x,0) = 0. 

4.5.5. Find the function u(x, t), defined for < x < tt and t > 0, which satisfies 
the following conditions: 

d 2 w d 2 u . . . , 

^ = ^2' «(0,t)=«(7T,t)=0, 

u(a;, 0) = 0, — (x, 0) = sin a; + sin 2a;. 

4.5.6. (For students with access to Mathematica) a. Find the first ten coeffi- 
cients of the Fourier sine series for 

h(x) = x — x 4 

by running the following Mathematica program 

f[n_] := 2 Nlntegrate [(x - xA4) Sin[n Pi x] , {x,0,l}]; 
b = Table [f [n] , {n,l,10}] 

b. Find the first ten terms of the solution to the initial value problem for a 
vibrating string, 

d 2 u d 2 u 

dt2=dx 2 > u(0,t)=«(7r,t)=0, 

u(x,0) = x — x 4 , ^— (x,0) = 0. 
at 

c. Construct a sequence of sketches of the positions of the vibrating string at 
the times ij = ih, where h — .1 by running the Mathematica program: 

vibstring = Table [ 
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Plot[ 

Sum[ b[n] Sin[n Pi x] Cos [n Pi t] , {n,l,10}], 
{x,0,l}, PlotRange -> {-1,1} 
], {t,0,l.,.l} 

] 

d. Select the cell containing vibstring and animate the sequence of graphics 
by running "Animate selected graphics," from the Cell menu. 



4.6 Heat flow in a circular wire 

The theory of Fourier series can also be used to solve the initial value problem 
for the heat equation in a circular wire of radius 1 which is insulated along the 
sides. In this case, we seek a function u(9,t), defined for 6 e R and t > such 
that 

1. u(6,t) satisfies the heat equation 

du 9 d 2 u 

1H= C W> (4 - 23) 

where c is a constant. 

2. u(0, t) is periodic of period 2tt in the variable 0; in other words, 

u(0 + 2ir,t)=u(6,t), for all and t. 

3. u(0,t) satisfies the initial condition u(0,O) = h(0), where h(0) is a given 
function, periodic of period 2n, the initial temperature of the wire. 

Once again the heat equation itself and the periodicity condition are homoge- 
neous and linear, so they must be dealt with first. Once we have found the 
general solution to the homogeneous conditions 

du 9 d 2 u , . 

-di = c W u (0 + 2 M = «(M), 

we will be able to find the particular solution which satisfies the initial condition 

u(0,O) = KB) 

by the theory of Fourier series. 

Thus we substitute u(0,t) = f(0)g(t) into the heat equation (4.23) to obtain 

f(0)g'(t) = c 2 f'{9)g{t) 

and separate variables: 

g'(t) _ c2 f"(0) 



g(t) m ' 
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A, or f"(6) = Xf(6). (4.25) 



The left-hand side of this equation does not depend on 9, while the right-hand 
side does not depend on t, so neither side can depend upon cither 9 or t, and 
we can write 

g'(t) = f"(9) = 
c 2 9(t) f{9) 

where A is a constant. We thus obtain two ordinary differential equations, 

J^pA, or 3 '(i) = Ac 2 5 (t), (4.24) 

and 

ng) 
m 

The periodicity condition u{9 + n,t) = u(9, t) becomes 

f(0 + 2w)g(t) = f(6)g(t), 
and if g(t) is not identically zero, we must have 

/(0 + 2tt) = /(0). 

Thus to find the nontrivial solutions to the homogeneous linear part of the 
problem requires us to find the nontrivial solutions to the problem: 

f"W = W {m) = Xm f(° + 2 ^ = f( 9 )- ( 4 - 26 ) 
We will call (4.11) the eigenvalue problem for the differential operator 

acting on the space V of functions which are periodic of period 2tt. 
As before, we need to consider three cases. 

Case 1: A = 0. In this case, the eigenvalue problem (4.26) becomes 

f"(9) = 0, f(9 + 2n) = f(9). 

The general solution to the differential equation is f(9) = a + b9, and 

a + b(9 + 2tt) = a + b{9) ^b = 0. 

Thus the only solution in this case is that where / is constant, and to be 
consistent with our Fourier series conventions, we write f(9) = a /2, where ao 
is a constant. 

Case 2: A > 0. In this case, the differential equation 

f"(6) = \f(6), or f"(9)-Xf(9)=0 
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has the general solution 
Note that 

a^O^ f{9) — > ±00 as 9 — > 00, 

while 

& 7^ /(0) -> ±00 as -» -00. 

Neither of these is consistent with the periodicity conditions f(9 + 2n) = f(6), 
so we conclude that a = b = 0, and we obtain no nontrivial solutions in this 
case. 

Case 3: A < 0. In this case, we set u = \/— A, and rewrite the eigenvalue 
problem as 

f"(9)+u; 2 f(9)=0, f(6 + 2w) = f(6). 

We recognize once again our old friend, the differential equation of simple har- 
monic motion, which has the general solution 

/(<?) = acos(Lu9) + bsin(Lu9) = Asin(uj(9 - 9 )). 

The periodicity condition f(9 + 2tt) = f(9) implies that u> = n, where n is an 
integer, which we can assume is positive, and 

f(9) = a n cos(n9) + b n sin(n9). 

Thus we see that the general solution to the eigenvalue problem (4.26) is 

A = and /(f) = y> 

or 

A = — n 2 where n is a positive integer and f(9) = a n cos(n9) + b n sm(n9). 
Now we need to find the corresponding solutions to (4.24) 

g'{t) = \c 2 g{t), 

for A = 0, — 1, —4, —9, . . . , — n 2 , .... As before, we find that the solution is 

g(t) = (constant)e~™ 2<;2 *, 

where c is a constant. Thus the product solutions to the homogeneous part of 
the problem are 

u (9, t) = Y> t) = K cos(n9) + b n sin(n0)]e _ " 2c2 *, 

where n = 1, 2, 3, ... . 
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Now we apply the superposition principle — an arbitrary superposition of 
these product solutions must again be a solution. Thus 

"(*>*) = y + CiK cos(nfl) + b n sin(n0)]e-" 2c2t (4-27) 

is a periodic solution of period 2tt to the heat equation (4.23). 

To finish the solution to our problem, we must impose the initial condition 

u(0,O) = h(6). 

But setting t = in (4.27) yields 

y + £~ =1 [a„cos(n0) + b n sin(n9)} = h{9), 

so the constants do, ai, . . . , £>i, . . . are just the Fourier coefficients of h{9). Thus 
the solution to our initial value problem is just (4.27) in which the constants a k 
and bk can be determined via the familiar formulae 

a k =- [ h{9)cosk9d6, b k = - [ h(9)sink9d9. 

^ J-TT ^ J-7T 

Note that as t — > oo the temperature in the circular wire approaches the constant 
value ao/2. 

Exercises: 

4.6.1. Find the function u(9, t), defined for < 9 < 2ir and t > 0, which satisfies 
the following conditions: 

— = — , u(6 + 2Tr,t) = u(6,t), u(9, 0) = 2 + sin 9 - cos 36*. 
You may assume that the nontrivial solutions to the eigenvalue problem 
f"(0) = Xf(0), f(9 + 27r) = f(9) 

are 

A = and f(9) = ^, 

and 

A = — n 2 and f{9) — a n cos n9 + b n sin n0, for n = 1, 2, 3, . . . , 
where the a„'s and &„'s are constants. 

4.6.2. Find the function u(9, t), defined for < 9 < 2ir and t > 0, which satisfies 
the following conditions: 

du d 2 u , n . , n . 
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u{9,0) = \9\, for 9 e [-tt.tt]. 



4.6.3. Find the function t), denned for < 6* < 27r and t > 0, which satisfies 
the following conditions: 

d 2 u d 2 u 

W = W «(* + 2 *> *) = «(*>*), 

u(0,O) = 2 + sin6>-cos30, — (0,O) = O. 



0t 



4.7 Sturm-Liouville Theory* 



We would like to be able to analyze heat flow in a bar even if the specific 
heat a{x), the density p{x) and the thermal conductivity k(x) vary from point 
to point. As we saw in Section 4.1, this leads to consideration of the partial 
differential equation 

'k{x)-\, (4.28) 



dt p(x)a(x) dx \ dx J ' 

where we make the standing assumption that p(x), <r(x) and k(x) are positive. 

We imagine that the bar is situated along the x-axis with its endpoints 
situated at x — a and x = b. As in the constant coefficient case, we expect that 
there should exist a unique function u(x,t), defined for a < x < b and t > 
such that 

1. u(x,t) satisfies the heat equation (4.28). 

2. u(x, t) satisfies the boundary condition u(a, t) = u(b, t) = 0. 

3. u(x,t) satisfies the initial condition u(x, 0) = h(x), where h(x) is a given 
function, defined for x G [a, b], the initial temperature of the bar. 

Just as before, we substitute u(x,t) — f(x)g(t) into (4.28) and obtain 

Once again, we separate variables, putting all the functions involving t on the 
left, all the functions involving x on the right: 

^- 1 d L(xAx^ 1 



g(t) p(x)a(x) dx \ dx J f(x) 

As usual, the two sides must equal a constant, which we denote by A, and our 
equation separates into two ordinary differential equations, 

g'(t) - Xg(t), (4.29) 
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and 



' ( K{x)^-{x))=\f{x). (4.80) 



p(x)a(x) dx \ dx 

Under the assumption that u is not identically zero, the boundary condition 
u(a, t) = u(b, t) = yields 

f(a) = f(b) = 0. 

Thus to find the nontrivial solutions to the homogeneous linear part of the 
problem, we need to find the nontrivial solutions to the boundary value problem: 

' '■•■:•<• )'? ■:•<•)) A/,:,-i. f(a) = = /(&). (4.31) 



p{x)a(x) dx \ dx 

We call this the eigenvalue problem or Sturm- Liouville problem for the differen- 
tial operator 

L 1 d ( K{x) d 

p(x)o~(x) dx \ dx 

which acts on the space Vo of well-behaved functions / : [a, b] — > R which 
vanish at the endpoints a and b. The eigenvalues of L are the constants A for 
which (4.31) has nontrivial solutions. Given an eigenvalue A, the corresponding 
eigenspace is 

W x = {/ e V Q : f satisfies (4.31)}. 

Nonzero elements of the eigenspaces are called eigenfunctions . 

If the functions p{x), o~(x) and k(x) are complicated, it may be impossible to 
solve this eigenvalue problem explicitly, and one may need to employ numerical 
methods to obtain approximate solutions. Nevertheless, it is reassuring to know 
that the theory is quite parallel to the constant coefficient case that we treated 
in previous sections. The following theorem, due to the nineteenth century 
mathematicians Sturm and Liouville, is proven in more advanced texts: 4 

Theorem. Suppose that p(x), cr(x) and k(x) are smooth functions which arc 
positive on the interval [a, b\. Then all of the eigenvalues of L are negative real 
numbers, and each eigenspace is one-dimensional. Moreover, the eigenvalues 
can be arranged in a sequence 

> Ai > A 2 > • • • > A„ > • • • , 

with X n — > — oo. Finally, every well-behaved function can be represented on 
[a, b] as a convergent sum of eigenfunctions. 

Suppose that fi(x), f2(x), ... , f n (x), . . . are eigenfunctions corresponding to 
the eigenvalues Ai, A2, ... , A„, .... Then the general solution to the heat 
equation (4.28) together with the boundary conditions u(a,t) = u(b,t) = is 



u(x,t) = ^2 Cnfn{x)e } ' 
n=0 



4 For further discussion, see Boyce and DiPrima, Elementary differential equations and 
boundary value problems, Seventh edition, Wiley, New York, 2001. 
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where the c„'s are arbitrary constants. 

To determine the c„'s in terms of the initial temperature h(x), we need a 
generalization of the theory of Fourier series. The key idea here is that the 
eigenspaces should be orthogonal with respect to an appropriate inner product. 
The inner product should be one which makes L like a symmetric matrix. To 
arrange this, the inner product that we need to use on Vq is the one defined by 
the formula 



g) = ( p{x)a{x)f{x)g{x)dx. 

J a 



Lemma. With respect to this inner product, eigenfunctions corresponding to 
distinct eigenvalues are perpendicular. 

The proof hinges on the fact that 

(L(/), ff ) = (/,L( 5 )}, bif,geV , 

so that if we thought of L as represented by a matrix, the matrix would be 
symmetric. This identity can be verified by integration by parts; indeed, 



(W),9)= I P (x)a(x)L(f)(x)g(x)dx = I ± ( K(xAx) ) g(x) 



b r b 

dx V dx 

b 



J a 



K{x)f x {x)f x {x)dx= ■■■ = {! M9)), 



where the steps represented by dots are just like the first steps, but in reverse 
order. 

It follows that if fi(x) and fj(x) are eigenfunctions corresponding to distinct 
eigenvalues A^ and Xj, then 

MfiJj) = (Mf),9) = (fMg)) = Wfufih 

and hence 

(\i-\ j )(f i ,f j ) = 0. 

Since Aj — Xj ^ 0, we conclude that fi and fj are perpendicular with respect to 
the inner product (•,•), as claimed. 

Thus to determine the c„'s, we can use exactly the same orthogonality tech- 
niques that we have used before. Namely, if we normalize the eigenfunctions 
so that they have unit length and are orthogonal to each other with respect to 
(,-,■>, then 

Cn — {h>: fn) i 

or equivalently, c„ is just the projection of h in the /„ direction. 



Example. We consider the operator 



d ( d 



dx \ dx 
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which acts on the space Vb of functions / : [1, e ir ] — > J? which vanish at the 
cndpoints of the interval [l,e 7r ]. To solve the eigenvalue problem, we need to 
find the nontrivial solutions to 

X Tx { X fx {x) ) = XI{X) ' = ° = /(e7r) - (4 - 32) 

We could find these nontrivial solutions by using the techniques we have learned 
for treating Cauchy-Euler equations. 

However, there is a simpler approach, based upon the technique of substitu- 
tion. Namely, we make a change of variables x — e z and note that since 

dx = e* dz , ' i£ and hence X ± = e«i£ = * 
ax e z az ax e z az dz 

Thus if we set 

f(z) = f(x) = f(e z ), 
the eigenvalue problem (4.32) becomes 

^(z) = Xf(z), /(0)=0 = /», 
a problem which we have already solved. The nontrivial solutions are 

A„ = — n 2 , fn(z) = sinnz, where n = 1,2, 

Thus the eigenvalues for our original problem are 

A n = — n 2 , for n = 1,2,..., 
and as corresponding cigenfunctions we can take 

f n (x) = sin(nloga;), 

where log denotes the natural or base e logarithm. The lemma implies that 
these cigenfunctions will be perpendicular with respect to the inner product 
(•,•), defined by 

(/,<?>= / -f(x)g(x)dx. (4.33) 
Ji x 

Exercises: 

4.7.1. Show by direct integration that if m ^ n, the functions 

fm{x) = sin(mlogx) and Jn{x) = sin(nlogx) 
are perpendicular with respect to the inner product defined by (4.33). 
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4.7.2. Find the function u(x,t), defined for 1 < x < and t > 0, which satisfies 
the following initial-value problem for a heat equation with variable coefficients: 

du d ( du\ ,„ , , _ . 

= «(c',*) = l 

u(x, 0) = 3 sin(log x) +7 sin(2 log x) — 2 sin(3 log x) . 

4. 7. 3. a. Find the solution to the eigenvalue problem for the operator 

which acts on the space Vo of functions / : [1, e n ] — > R which vanish at the 
cndpoints of the interval [1, e^]. 

b. Find the function u(x,t), defined for 1 < x < e" and t > 0, which satisfies 
the following initial-value problem for a heat equation with variable coefficients: 

m =x dx{ x dx)-^ u ^ = < e '*) = ' 

u(x, 0) = 3 sin(log x) +7 sin(2 log x) — 2 sin(3 log x) . 



4.8 Numerical solutions to the eigenvalue prob- 
lem* 

We can also apply Sturm-Liouville theory to study the motion of a string of 
variable mass density. We can imagine a violin string stretched out along the 
x-axis with endpoints at x — and x — 1 covered with a layer of varnish which 
causes its mass density to vary from point to point. We could let 

p(x) = the mass density of the string at x for < x < 1. 

If the string is under constant tension T, its motion might be governed by the 
partial differential equation 

d 2 u T d 2 u 

(4.34) 



dt 2 p(x) dx 2 ' 

which would be subject to the Dirichlct boundary conditions 

u(0,t) = = u(l,t), forallt>0. (4.35) 

It is natural to try to find the general solution to (4.34) and (4.35) by sepa- 
ration of variables, letting u(x, t) = f(x)g(t) as usual. Substituting into (4.34) 
yields 

t< \ "m T ft ^ m 9"(t) T f'(x) 

f(-)9(t) = W) f(x) 9 (t), or ^-^y-^- 
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The two sides must equal a constant, denoted by A, and the partial differential 
equation separates into two ordinary differential equations, 

^- )f "{ x ) = Xf{x), g"(t) = \g(t). 

The Dirichlct boundary conditions (4.35) yield /(0) = = /(l). Thus / must 
satisfy the eigenvalue problem 

L(/) = A/, /(0)=0 = /(l), where L= ''' 



p{x) dx 2 

Although the names of the functions appearing in L are a little different than 
those used in the previous section, the same Theorem applies. Thus the eigen- 
values of L are negative real numbers and each eigenspace is one-dimensional. 
Moreover, the eigenvalues can be arranged in a sequence 

> Ai > A 2 > • • • > A„ > • • • , 

with A„ — > — oo. Finally, every well-behaved function can be represented on 
[a, b] as a convergent sum of eigcnfunctions. If /i(x), f2( x ), ■ ■ ■ , fn( x ), ■ ■ ■ are 
cigcnfunctions corresponding to the eigenvalues Ai, A 2 , ... , A„, ... . Then the 
general solution to (4.34) and (4.35) is 

oo 

u(x, t) = f n (x) a n cos(v/-A„t) + b n sm(y/-X n t) 

n=0 

where the a„'s and 6„'s are arbitrary constants. Each term in this sum represents 
one of the modes of oscillation of the vibrating string. 

In constrast to the case of constant density, it is usually not possible to find 
simple explicit eigenfunctions when the density varies. It is therefore usually 
necessary to use numerical methods. 

The simplest numerical method is the one outlined in §4.3. For < i < n, 
we let Xi — i/n and 

Ui(t) = u(xi,t) = the displacement at a;, at time t. 

Since Uo(t) = = u n (t) by the boundary conditions, the displacement at time 
t is approximated by 

/ Ul (t) \ 
u(t)= U2{t) 

\U n -i(t)J 

a vector-valued function of one variable. The partial derivative 

d 2 u . . d 2 u 

-^p 1S approximated by 
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Figure 4.1: Shape of the lowest mode when p = l/(x + .1). 



and as we saw in §4.3, the partial derivative 
8 2 u 



dx 2 



is approximated by n Pu, 



where 



/ -2 1 
1 -2 1 
1 -2 



\ 




\ • ••• -2 ) 
Finally, the coefficient (T/p(x)) can be represented by the diagonal matrix 

\ 



Q 



( T/p( Xl ) 
T/p(x 2 ) 
T/p(x 3 ) 



V 











T/p(x n ^) J 



Putting all this together, we find that our wave equation with variable mass den- 
sity is approximated by a second order homogeneous linear system of ordinary 
differential equations 



dt 2 



= Au, where A = n 2 QP. 



The eigenvalues of low absolute value are approximated by the eigenvalues of A, 
while the eigenfunctions representing the lowest frequency modes of oscillation 
are approximated eigenvectors corresponding to the lowest eigenvalues of A. 

For example, we could ask the question: What is the shape of the lowest 
mode of oscillation in the case where p(x) = l/(x+ .1)? To answer this question, 
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we could utilize the following Mathematica program (which is quite similar to 
the one presented in Exercise 2.5.2): 



n := 100; rho [x_] := l/(x + .1); 

m := Table [Max [2-Abs [i-j] ,0] , { i,n-l } ,{ j,n-l } ]; 
p := m - 4 IdentityMatrix [n-1] ; 

q := DiagonalMatrix [Table [( 1/rho [i/n] ) , { i,l,n-l } ]]; 
a := nA2 q.p; eigenvec = Eigenvectors [N [a] ] ; 
ListPlot [eigenvec [ [n-1] ] ] 

If we run this program we obtain a graph of the shape of the lowest mode, 
as shown in Figure 4.1. Note that instead of approximating a sine curve, our 
numerical approximation to the lowest mode tilts somewhat to the left. 
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Chapter 5 



PDE's in Higher 
Dimensions 



5.1 The three most important linear partial dif- 
ferential equations 

In higher dimensions, the three most important linear partial differential equa- 
tions are Laplace's equation 

d 2 u d 2 u d 2 u 
dx 2 dy 2 dz 2 

the heat equation 

du 
dt 

and the wave equation, 

d 2 u 2 ( d 2 u d 2 u d 2 u\ 
dt 2 \ dx 2 dy 2 dz 2 J ' 

where c is a nonzero constant. Techniques developed for studying these equa- 
tions can often be applied to closely related equations. 

Each of these three equations is homogeneous linear, that is each term con- 
tains u or one of its derivatives to the first power. This ensures that the principle 
of superposition will hold, 

U\ and u 2 solutions C\U\ + c 2 u 2 is a solution, 

for any choice of constants c\ and c 2 . The principle of superposition is essential 
if we want to apply separation of variables and Fourier analysis techniques. 

In the first few sections of this chapter, we will derive these partial differen- 
tial equations in several physical contexts. We will begin by using the divergence 



/ d 2 u d 2 u d 2 u\ 
\ dx 2 dy 2 dz 2 J ' 



115 



theorem to derive the heat equation, which in turn reduces in the steady-state 
case to Laplace's equation. We then present two derivations of the wave equa- 
tion, one for vibrating membranes and one for sound waves. Exactly the same 
wave equation also describes electromagnetic waves, gravitational waves, or wa- 
ter waves in a linear approximation. It is remarkable that the principles devel- 
oped to solve the three basic linear partial differential equations can be applied 
in so many contexts. 

In a few cases, it is possible to find explicit solutions to these partial dif- 
ferential equations under the simplest boundary conditions. For example, the 
general solution to the one-dimensional wave equation 

d 2 u _ 

for the vibrations of an infinitely long string, is 

u(x, t) = f(x + ct) + g(x — ct) , 

where / and g are arbitrary well-behaved functions of a single variable. 

Slightly more complicated cases require the technique of "separation of vari- 
ables" together with Fourier analysis, as we studied before. Separation of vari- 
ables reduces these partial differential equations to linear ordinary differential 
equations, often with variable coefficients. For example, to find the explicit so- 
lution to the heat equation in a circular room, we will see that it is necessary 
to solve Bessel's equation. 

The most complicated cases cannot be solved by separation of variables, 
and one must resort to numerical methods, together with sufficient theory to 
understand the qualitative behaviour of the solutions. 

In the rest of this section, we consider our first example, the equation gov- 
erning heat flow through a region of (x, y, z)-space, under the assumption that 
no heat is being created or destroyed within the region. Let 

u(x, y, z, t) = temperature at (a;, y, z) at time t. 

If a(x,y,z) is the specific heat at the point (x,y,z) and p(x,y,z) is the 
density of the medium at (x,y,z), then the heat within a given region D in 
(x, y, z)-space is given by the formula 

Heat within D = J J J p(x,y, z)cr(x,y, z)u(x,y, z,t)dxdydz. 

If no heat is being created or destroyed within D, then by conservation of energy, 
the rate at which heat leaves D equals minus the rate of change of heat within 
D, which is 



d_ 
~~dt 



paudxdydz 

D 



pa—dxdydz, (5-1) 

d 



by differentiating under the integral sign. 
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On the other hand, heat flow can be represented by a vector field F(x, y, z, t) 
which points in the direction of greatest decrease of temperature, 

F(x,y,z,t) = -n(x,y,z)(Vu)(x,y,z,t), 

where n(x,y,z) is the so-called thermal conductivity of the medium at (x,y, z). 
Thus the rate at which heat leaves the region D through a small region in its 
boundary of area dA is 

-(kVu) • NdA, 

where N is the unit normal which points out of D. The total rate at which heat 
leaves D is given by the flux integral 



- / / («V«) • NdA, 



where dD is the surface bounding D. It follows from the divergence theorem 
that 

Rate at which heat leaves D = — j j j V • {nVu)dxdydz. (5.2) 

From formulae (5.1) and (5.2), we conclude that 

J J J P<y-^dxdydz — J J J V • (nVu)dxdydz. 

This equation is true for all choices of the region D, so the integrands on the 
two sides must be equal: 

p(x,y,z)a(x,y,z) — (x,y,z,t) = V • ( K Wu)(x,y, z,t). 



Thus we finally obtain the heat equation 
du 1 



dt p(x,y,z)a(x,y,z) 



V -(K(x,y,z)(Vu)) 



In the special case where the region D is homogeneous, i.e. its properties are 
the same at every point, p(x,y,z), a(x,y,z) and n{x,y,z) are constants, and 
the heat equation becomes 



du K 
dt pa 



d 2 u d 2 u d 2 u 
dx 2 dy 2 dz 2 



If we wait long enough, so that the temperature is no longer changing, the 
"steady-state" temperature u(x,y,z) must satisfy Laplace's equation 

d 2 u d 2 u d 2 u 
dx 2 dy 2 dz 2 
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If the temperature is independent of z, the function u(x,y) — u(x,y,z) must 
satisfy the two-dimensional Laplace equation 



d 2 u ^ d 2 u ^ 
dx 2 dy 2 



Exercises: 



5.1.1. For which of the following differential equations is it true that the super- 
position principle holds? 

d 2 u d 2 u d 2 u 
dx 2 dy 2 dz 2 

d 2 u d 2 u d 2 u 
ox z oy z oz z 

d 2 u d 2 u d 2 u 2 
dx 2 dy 2 dz 2 

d 2 u d 2 u d 2 u _ x 
dx 2 dy 2 dz 2 

d 2 u d 2 u d 2 u x 
dx 2 dy 2 dz 2 

Explain your answers. 

5.1.2. Suppose that a chemical reaction creates heat at the rate 

X(x, y, z)u(x, y, z, t) + v{x, y, z) , 
per unit volume. Show that in this case the equation governing heat flow is 
Ou 

p(x,y,z)a(x,y,z)— = X(x,y, z)u(x,y, z,t) + v{x,y,z) + V • (k(x, y, z)(Vn)) . 

5.2 The Dirichlet problem 

The reader will recall that the space of solutions to a homogeneous linear second 
order ordinary differential equation, such as 

d 2 u du 

— +P^ + 9 u = 

is two-dimensional, a particular solution being determined by two constants. By 
contrast, the space of solutions to Laplace's partial differential equation 

d 2 u d 2 u n 

dx~ 2+ W = ° ^ 
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is m/mi£e-dimensional. For example, the function u(x, y) = x 3 — 3xy 2 is a 
solution to Laplace's equation, because 

du 2 2 9 2 u 



<9u „ <9 2 u 
— = -Qxy, — 
dy dy 



Gxy, = -6x, 



and hence 



d 2 u d 2 u 

— + — = 6x - 6x = 0. 
ax z ay" 2 



Similarly, 

u(x,y)=7, u(x, y) = x 4 — 6x 2 y 2 + y A , and u(x, y) — e x sin y 

are solutions to Laplace's equation. A solution to Laplace's equation is called 
a harmonic function. It is not difficult to construct infinitely many linearly 
independent harmonic functions of two variables. 

Thus to pick out a particular solution to Laplace's equation, we need bound- 
ary conditions which will impose infinitely many constraints. To see what 
boundary conditions are natural to impose, we need to think of a physical 
problem which leads to Laplace's equation. Suppose that u(x,y) represents 
the steady-state distribution of temperature throughout a uniform slab in the 
shape of a region D in the (x, y)-planc. If we specify the temperature on the 
boundary of the region, say by setting up heaters and refrigerators controlled by 
thermostats along the boundary, we might expect that the temperature inside 
the room would be uniquely determined. We need infinitely many heaters and 
refrigerators because there are infinitely many points on the boundary. Spec- 
ifying the temperature at each point of the boundary imposes infinitely many 
constraints on the harmonic function which realizes the steady-state tempera- 
ture within D. 

The Dirichlet Problem for Laplace's Equation. Let D be a bounded 
region in the (x, y)-plane which is bounded by a curve dD, and let <j> : dD — > R 
be a continuous function. Find a harmonic function u : D — > R such that 

u{x, y) = 4>{x, y), for (x, y) € dD. 



Our physical intuition suggests that the Dirichlet problem will always have a 
unique solution. This is proven mathematically for many choices of boundary 
in more advanced texts on complex variables and partial differential equations. 
The mathematical proof that a unique solution exists provides evidence that 
the mathematical model we have constructed for heat flow may in fact be valid. 

Our goal here is to find the explicit solutions in the case where the region D 
is sufficiently simple. Suppose, for example, that 

D = {{x, y) £ R 2 : < x < a, < y < b}. 
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Suppose, moreover that the function 4> ■ dD — > R vanishes on three sides of 3D, 
so that 

<t>{0,y) = cj>{a,y) = <f>{x,0)=0, 

while 

4>{x,b) = f(x), 

where f{x) is a given continuous function which vanishes when x = and x — a. 

In this case, we seek a function u(x, y), defined for < x < a and < y < b, 
such that 

1. u(x,y) satisfies Laplace's equation 

d 2 u d 2 u „ . ,, 

dx-* + W=°- ^ 

2. u(x,y) satisfies the homogeneous boundary conditions u(Q,y) — u(a,y) = 
u(x,0) = 0. 

3. u(x,y) satisfies the nonhomogeneous boundary condition u(x,b) = f(x), 
where f(x) is a given function. 

Note that the Laplace equation itself and the homogeneous boundary con- 
ditions satisfy the superposition principle — this means that if u\ and satisfy 
these conditions, so does c\U\ + C2U2, for any choice of constants c\ and C2. 

Our method for solving the Dirichlet problem consists of two steps: 

Step I. We find all of the solutions to Laplace's equation together with the 
homogeneous boundary conditions which are of the special form 

u(x,y)=X(x)Y(y). 

By the superposition principle, an arbitrary linear superposition of these solu- 
tions will still be a solution. 

Step II. We find the particular solution which satisfies the nonhomogeneous 
boundary condition by Fourier analysis. 

To carry out Step I, we substitute u(x, y) = X(x)Y(y) into Laplace's equation 
(5.3) and obtain 

X"(x)Y(y) + X(x)Y"(y) = 0. 

Next we separate variables, putting all the functions involving x on the left, all 
the functions involving y on the right: 

X"(x) Y"{y) 
X(x) ~ Y(y) ■ 

The left-hand side of this equation does not depend on y, while the right-hand 
side does not depend on x. Hence neither side can depend upon cither x or y. 
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In other words, the two sides must equal a constant, which we denote by A, and 
call the separating constant, as before. Our equation now becomes 

X"(x) _ Y"(y) = 
X{x) Y{y) 

which separates into two ordinary differential equations, 

X"(x) = XX(x), (5.5) 

and 

Y"(y) = -\Y(y). (5.6) 
The homogeneous boundary condition u(0,y) — u(a,y) = imply that 
X(0)Y(y) = X(a)Y(y) = 0. 
If Y(y) is not identically zero, 

X(0) = X(a) = 0. 

Thus we need to find the nontrivial solutions to a boundary value problem for 
an ordinary differential equation: 

X"(x) = ^{X{x)) = XX(x), X(0) = = X(a), 
which we recognize as the eigenvalue problem for the differential operator 




acting on the space V of functions which vanish at and a. We have seen before 
that the only nontrivial solutions to equation (5.5) are constant multiples of 

X(x) = sm(mrx/a), with A = — (mr/a) 2 , n= 1,2,3,.... 

For each of these solutions, we need to find a corresponding Y{y) solving 
equation (5.6), 

Y"(y) = -XY(y), 

where A = —(nir/a) 2 , together with the boundary condition Y(0) = 0. The 
differential equation has the general solution 

Y{y) = Acosh(mry /a) + B smh(mry / a) , 

where A and B are constants of integration, and the boundary condition Y(0) = 
implies that A = 0. Thus we find that the nontrivial product solutions to 
Laplace's equation together with the homogeneous boundary conditions are con- 
stant multiples of 

u n (x, y) = sm{mtx/a) sinh(n7ry/a). 
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The general solution to Laplace's equation with these boundary conditions is a 
general superposition of these product solutions: 

u(x, y) = Bi sin(7ra/a) smh(ny/a) + B 2 sm(2nx/a) sinh(27rj//a) + . . . . (5.7) 

To carry out Step II, we need to determine the constants Bi,B 2 ,... which 
occur in (5.7) so that 

u(x,b) = f(x). 
Substitution of y = b into (5.7) yields 

f(x) = Bi sm(nx/a) sinh(7r6/a) + B 2 sin(27rx/a) sinh(27r6/a)+ 

. . . + B k sin(27rfc/o) sinh(A;7r6/a) + . . . . 

We see that Bk smh(knb/a) is the fc-th coefficient in the Fourier sine scries for 
/(*)• 

For example, if a = b = n, and f(x) = 3 sin x + 7sin2x, then we must have 
f(x) = Bi sin(x) sinh(7r) + B 2 sin(2x) sinh(27r), 

and hence 

Si = ^--, B 2 = , J B k ={) for fe = 3,.... 

smh(7r) smh(27r) 

Thus the solution to Dirichlet's problem in this case is 

3 7 

u(x, y) = . - , , sin x sinh y + . - — . sin 2x sinh 2y. 
smh(7r) smh(27r) 

Exercises: 

5.2.1. Which of the following functions are harmonic? 
a - f{x,y) = x 2 + y 2 . 

b. f(x,y) =x 2 -y 2 . 
c f{x,y) = e x cosy. 
d. f(x,y) =x 3 -3xy 2 . 

5.2.2. a. Solve the following Dirichlet problem for Laplace's equation in a square 
region: Find u(x, y),0<x<Tr,0<y<n, such that 

3 2 u d 2 u 

^2+^2=°' u(0,y)=u(7r,») = 0, 
u(x, 0) = 0, u(x, 7r) = sin x — 2 sin 2x + 3 sin 3x. 
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Figure 5.1: Graph of u(x, y) — sin ^^ sin x sinh y + sinll ( 271 -) sin 2x sinh 2y. 



b. Solve the following Dirichlet problem for Laplace's equation in the same 
square region: Find u(x, y), < x < tt, < y < tt, such that 



+ 1^2-0, u(0,y) = 0, 



dx 2 dy 2 

u(ir, y) = sin 2y + 3 sin4y, u(x, 0) = = u(x, tt). 



c. By adding the solutions to parts a and c together, find the solution to the 
Dirichlet problem: Find u(x, y), < x < tt, < y < tt, such that 

d 2 u d 2 u 

dx~ 2+ W = Q ' W( °' 2/) = ' 
u(tt, y) — sin 2j/+3 sin Ay, u(x, 0) = 0, u(x, tt) = sin x— 2 sin 2x+3 sin 3x. 

5.3 Initial value problems for heat equations 

The physical interpretation behind the heat equation suggests that the following 
initial value problem should have a unique solution: 

Let D be a bounded region in the {x, ?/)-plane which is bounded by a piece- 
wise smooth curve 3D, and let h : D — > R be a continuous function which 
vanishes on dD. Find a function u(x,y,t) such that 
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1. u satisfies the heat equation 

du of d 2 u d 2 u\ 

m= c [d^ + WJ- (5 ' 8) 

2. u satisfies the "Dirichlet boundary condition" u(x,y,t) = 0, for (x,y) e 
3D. 

3. u satisfies the initial condition u(x, y, 0) = h(x, y). 

The first two of these conditions are homogeneous and linear, so our strategy is 
to treat them first by separation of variables, and then use Fourier analysis to 
satisfy the last condition. 

In this section, we consider the special case where 

D = {(x, y) e R 2 : < x < a, < y < 6}, 

so that the Dirichlet boundary condition becomes 

u(0, y, t) = u(a, y, t) = u(x, 0, t) = u(x, b, t) = 0. 

In this case, separation of variables is done in two stages. First, we write 

u(x,y,t) = f(x,y)g(t), 

and substitute into (5.8) to obtain 

/(,, y ),'(t)^ 2 (g + g)^). 

Then we divide by c 2 f(x,y)g(t), 



c 2 g(t) f{x,y) \dx 2 dy 

The left-hand side of this equation does not depend on x or y while the right 
hand side does not depend on t. Hence neither side can depend on x, y, or t, so 
both sides must be constant. If we let A denote the constant, we obtain 

1 VM = a 1 ■ a2f 



c 2 g(t)' f(x,y) \dx 2 dy 2 , 

This separates into an ordinary differential equation 

g'(t) = c 2 Xg(t), (5.9) 
and a partial differential equation 



2 f «2 . , 

A/, (5.10) 



\ dx 2 dy 2 
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called the Helmholtz equation. The Dirichlct boundary condition yields 

f(0,y) = f(a,y) = f(x,0) = f(x,b)=0. 

In the second stage of separation of variables, we set f(x,y) = X(x)Y(y) 
and substitute into (5.10) to obtain 

X"(x)Y(y) +X(x)Y"(y) = XX(x)Y(y), 

which yields 

X"(x) Y"(y) 
X(x) + Y(y) 

or 

X"(x) Y"(y) 
X(x) Y(y) ■ 

The left-hand side depends only on x, while the right-hand side depends only 
on y, so both sides must be constant, 

x"(x) \ _ 

X(x) " Y(y) ■ 

Hence the Helmholtz equation divides into two ordinary differential equations 

X"(x) = fiX(x), Y"(y) = vY(y), where fi + v = A. 
The "Dirichlet boundary conditions" now become conditions on X(x) and Y(y): 

X(0) = X(a) = 0, y(0) = Y(b) = 0. 
The only nontrivial solutions are 

X(x) = sin(m7rx/a), with /i = — (unr/a) 2 , m= 1,2,3,..., 

and 

Y(y) = sm(mry/b), with v=—(nn/b) 2 , n= 1,2,3,.... 
The corresponding solutions of the Helmholtz equation are 

fmn(x,y) = sin(irnvx/a)sm(niry/b), with X mn = — (ran /a) 2 — (mr/b) 2 . 

For any given choice of m and n, the corresponding solution to (5.9) is 

g(t) = e~ c2Xmnt = e - c2 (( m7r / a ) 2 +(™ 7r / b ) 2 )*. 

Hence for each choice of m and n, we obtain a product solution to the heat 
equation with Dirichlct boundary condition: 

u m , n {x,y,t) = sin(m7ra ; / a )sin(n7ry/6)e~ c2((m7r/Q)2+( ™ 7r/h)2) *. 
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The general solution to the heat equation with Dirichlet boundary conditions is 
an arbitrary superposition of these product solutions, 

oo 

u{x,y,t)= b mn sm(mnx/a)sm(nTTy/b)e- c2( - ( - m7r/a '> 2+( > njr/b)2)t . (5.11) 

m,n— 1 

To find the constants b mn appearing in (5.11), we need to apply the initial 
condition u(x,y,0) = h(x,y). The result is 

oo 

h(x,y) = b mn sm(mirx/a) sm(niry/b), (5-12) 

m,n— 1 

expressing the fact that the b mn 's are the coefficients of what is called the double 
Fourier sine series of h. As in the case of ordinary Fourier sine series, the b mn 's 
can be determined by an explicit formula. To determine it, we multiply both 
sides of (5.12) by (2/a) sin(p7ra/a), where p is a positive integer, and integrate 
with respect to x to obtain 



2 f a 

- / h(x,y) 
a Jo 



sm(pTTx/a)dx 



r 2 



= E * 



m,n—l 



a 



sm(pirx/a) sm(mTTx/a)dx 



ii 



sin(mry/b). 



The expression within brackets is one if p — m and otherwise zero, so 



/ h(x, y) sin(pirx/a)dx — \] b pn sm(mry/b). 



Next we multiply by (2/6) sm(qTry/b), where q is a positive integer, and integrate 
with respect to y to obtain 



2 2 r " rb 
ab ./,, ./„ 



h{x, y) sm{pitx/a) shi(qiry/b)dxdy 



n=l 



2 ' 6 



sin(n7n//&) sin(qiry /b)dy 



The expression within brackets is one if g = n and otherwise zero, so we finally 
obtain 



2 2 r a 

b P q—--r / / h(x,y)sin(pirx/a)sin(qiry/b)dxdy. (5.13) 
a ° Jo Jo 

Suppose, for example, that c = 1, a = 6 = 7r and 

/i(a:, y) = sin x sin y + 3 sin 2x sin j/ + 7 sin 3x sin 2y. 
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In this case, we do not need to carry out the integration indicated in (5.13) 
because comparison with (5.12) shows that 

611 = 1, 621 = 3, 632 = 7, 

and all the other b mn 's must be zero. Thus the solution to the initial value 
problem in this case is 

u(x, y, t) = sin a; sin ye -2 * + 3 sin 2x sin ye~ 5t + 7sin3ccsin2ye -13 *. 

Here is another example. Suppose that a = b = tt and 

h(x,y) =p(x)q(y), 

where 



p(x) = I 
In this case, 



x, for < x < 7r/2, 

TT — X, for 7r/2 < X < TT, 



q(y) 



y, for < y < tt/2, 

tt — y, for tt/2 <y <tt. 



b mn = 0^ J J p(x)q(y) sin mxdx 

p(x) sin mxdx 

= I — I / x sin mxdx + (tt — x) sin mxdx 
V 71 / J a Jtt/2 



sin nydy 
q(y) srnnydy 



/ y sin nydy + (tt — y) sin nydy 

JO Jtt/2 



The integration can be carried out just like we did in Section 3.3 to yield 
2 



m" 



■ sin(m7r/2) 



■ sin(n7r/2) 



16 

- ^2 

Thus in this case, we see that 



— ~ — ^ sin(m7r/2) sin(n7r/2) 
m n z 



16 



m,n— 1 



~tt \ sin(m7r/2) sin(n7r/2) 
m z n z 



sin mi sin my, 



and hence 



. 16 \ 
u(x,y,t) = -j 2^ 



m,n— 1 



— - — sin(m7r/2) sin(n7r/2) 
m n 1 



sin mi sin mye 
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Exercises: 

5.3.1. Solve the following initial value problem for the heat equation in a square 
region: Find u(x, y, t), where 0<x<7r, 0<y<7r and t > 0, such that 

du d 2 u d 2 u 
dt dx 2 dy 2 ' 

u{x, 0, t) = u(x, 7r, t) = u(0, y, t) — u(tt, y, t) = 0, 
u(x, y, 0) = 2 sin x sin y + 5 sin 2x sin y. 
You may assume that the nontrivial solutions to the eigenvalue problem 

0(x,y) + 0(z,tf) - A/foi,), /(x.O) = /(s.tt) = /(0,y) = /(tt,i/) = 0, 
are of the form 

A = —to 2 — n 2 , f{x,y) = b mn smmxsmny, 
for to = 1, 2, 3, . . . and n = 1,2,3,..., where b mn is a constant. 

5.3.2. Solve the following initial value problem for the heat equation in a square 
region: Find u(x, y, t), where 0<x<7r, 0<y<7r and t > 0, such that 

du d 2 u d 2 u 
dt dx 2 dy 2 ' 

u(x, 0, t) = u(x, 7r, t) — u(0, y, t) = u(ir, y, t) = 0, 
u(x,y,0) = 2(smx)y(n - y). 

5.3.3. Solve the following initial value problem in a square region: Find u(x, y, t), 
where 0<a;<7r, 0<y<7r and t > 0, such that 

1 du _ d 2 u d 2 u 
2t+l dt ~ dx 2 + dy 2 ' 

u(x, 0, t) = u(x, 7r, t) = u(0, y, t) — u(ir, y, t) = 0, 
u(x, y, 0) = 2 sin x sin y + 3 sin 2x sin y. 

5.3.4. Find the general solution to the heat equation 

du d 2 u d 2 u 
dt dx 2 dy 2 

subject to the boundary conditions 

u(x, 0, t) — u(0, y, t) = u(ir, y, t) = 0, 
u(x, Tr,t) = sin x — 2 sin 2x + 3 sin 3x. 
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5.4 Two derivations of the wave equation 



We have already seen how the one-dimensional wave equation describes the 
motion of a vibrating string. In this section we will show that the motion of a 
vibrating membrane is described by the two-dimensional wave equation, while 
sound waves are described by a three-dimensional wave equation. In fact, we 
will see that sound waves arise from "linearization" of the nonlinear equations 
of fluid mechanics. 

The vibrating membrane. Suppose that a homogeneous membrane is fas- 
tened down along the boundary of a region D in the (x, y)-plane. Suppose, 
moreover, that a point on the membrane can move only in the vertical direc- 
tion, and let u{x,y,t) denote the height of the point with coordinates (x,y) at 
time t. 

If p denotes the density of the membrane (assumed to be constant), then 
by Newton's second law, the force F acting on a small rectangular piece of 
the membrane located at (x, y) with sides of length dx and dy is given by the 
expression 

F = p-^-(x,y)dxdyk. 

Suppose the force displaces the membrane from a given position u(x, y) to a 
new position 

u(x,y) +rj(x,y), 

where r](x, y) and its derivatives are very small. Then the total work performed 
by the force F will be 

d it 

F • r)(x, y)k = tj(x, y)p-^j(x, y)dxdy. 

Integrating over the membrane yields an expression for the total work performed 
when the membrane moves through the displacement 77: 

f f d 2 u 
Work = J J i](x,y)p—(x,y)dxdy. (5.14) 

On the other hand, the potential energy stored in the membrane is propor- 
tional to the extent to which the membrane is stretched. Just as in the case of 
the vibrating string, this stretching is approximated by the integral 



T 

Potential energy = — 



D 



du\ 2 ( du N 



dx ) \dy 



dxdy, 



where T is a constant, called the tension in the membrane. Replacing u by u + n 
in this integral yields 



T 

New potential energy = — 



D 



du\ 2 ( du x 



dx ) \dy 



dxdy 
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+T 



D I 



dx J \dx J \dy J \dy 



dxdy 



T 
+ 2 



drj\ 2 ( dt]^ 



dx J \dy 



dxdy. 



If we neglect the last term in this expression (which is justified if r\ and its 
derivatives are assumed to be small), we find that 

New potential energy — Old potential energy = J J TVu • \7rjdxdy. 

It follows from the divergence theorem in the plane and the fact that r] vanishes 
on the boundary dD that 

/ / TVu ■ Vrjdxdy +11 TrjV ■ Vudxdy = / TrjVu ■ Nds = 0, 

J JD J JD JdD 

and hence 

Change in potential — — I I r)(x,y)T(V ■ Vu)(x,y)dxdy. (5.15) 



The work performed must be minus the change in potential energy, so it 
follows from (5.14) and (5.15) that 

/ J D v ( x,y } p ^( x,y ^ dxdy = J j D V(x,y)T(\7 -Vu)(x,y)dxdy. 
Since this equation holds for all choices of 77, it follows that 



which simplifies to 



or equivalcntly 



^ = rv.v«, 



d 2 u 2 „ „ , 2 T 

— — = c V • Vu, where c = — , 
at 1 p 



d 2 u 2 / d 2 u d 2 u 
dt 2 \ dx 2 dy 2 
This is just the wave equation. 

The equations of a perfect gas and sound waves. Next, we will describe 
Euler's equations for a perfect gas in (x, y, ^-space. 1 Euler's equations are 
expressed in terms of the quantities, 

v(x, y, z, t) = (velocity of the gas at (x, y, z) at time t), 



1 For a complete derivation of these equations, see Chapter 9 of Alexander Fetter and John 
Walecka, Theoretical mechanics of particles and continua, McGraw-Hill, New York, 1980, or 
Chapter 1 of Alexandre Chorin and Jerrold Marsden, A mathematical introduction to fluid 
mechanics, third edition, Springer, 1993. 



130 



p(x,y,z,t) = (density at (x,y,z) at time t), 
p(x,y,z,t) — (pressure at (x,y, z) at time t). 
The first of the Euler equations is the equation of continuity, 

^+V-(pv) = 0. (5.16) 

To derive this equation, we represent the fluid flow by the vector field 

F = pv, 

so that the surface integral 

IL FNdA 

represents the rate at which the fluid is flowing across S in the direction of N. 
We assume that no fluid is being created or destroyed. Then the rate of change 
of the mass of fluid within D is given by two expressions, 

J J J -J^( x -,y, z -,t)dxdydz and J J F ' Nc ^' 

which must be equal. It follows from the divergence theorem that the second of 
these expressions is 



and hence 



ill ^ '^^ x, y ,z ^^ x< ^y^ z ' 



^-dxdyd: = - / / / V • Fd.rdt/dz 



ID UL J J JD 

Since this equation must hold for every region D in (x, y, z)-space, we conclude 
that the integrands must be equal, 

! = -V.F = -V.(pv), 

which is just the equation of continuity. 

The second of the Euler equations is simply Newton's second law of motion, 

(mass density) (acceleration) = (force density). 

We make an assumption that the only force acting on a fluid element is due to 
the pressure, an assumption which is not unreasonable in the case of a perfect 
gas. In this case, it turns out that the pressure is defined in such a way that 
the force acting on a fluid element is minus the gradient of pressure: 

Forces -Vp(x(t), y(t), z(t), t). (5.17) 
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The familiar formula Force = Mass x Acceleration then yields 

P j t (v(x(t),y(t),z(t),t)) = -Wp{x{t),y{t),z{t),t). 



Using the chain rule. 



dv dv dx dv dy dv dz dv dt 
dt dx dt dy dt dz dt dt dt'' 

we can rewrite this equation as 

| + (vV)v=-V (5.18) 

Note that this equation is nonlinear because of the term (v • V)v. 

To finish the Euler equations, we need an equation of state, which relates 
pressure and density. The equation of state could be determined by experiment, 
the simplest equation of state being 

p=aV, (5.19) 

where a 2 and 7 are constants. (An ideal monatomic gas has this equation of 
state with 7 = 5/3.) 

The Euler equations (5.16), (5.18), and (5.19) arc nonlinear, and hence quite 
difficult to solve. However, one explicit solution is the case where the fluid is 
motionless, 

P = Pa, P = Pa, v = 0, 

where po and po satisfy 

2 7 

Pa = a p(j. 

Linearizing Euler's equations near this explicit solution gives rise to the linear 
differential equation which governs propagation of sound waves. 
Let us write 

P = Pa + P',P = Po + P', v = v', 

where p' , p' and v' are so small that their squares can be ignored. 
Substitution into Euler's equations yields 

% + PoV- (v') = 0, 

-qT = Vp', 

dt po 

and 

p'^a^po^-V = cV, 
where c 2 is a new constant. It follows from these three equations that 

d2p ' - -p fv.^Uv.Vp-c 2 V.Vp'. 



dt 2 ru v dt 
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Thus p' must satisfy the three-dimensional wave equation 




c 2 V • Vp' 



2 



d 2 P ' d 2 P ' d 2 P ' 



) 



(5.20) 



dt 2 



= c 



dx 2 dy 2 dz 2 



If the sound wave p' is independent of z, (5.20) reduces to 




d 2 p' d 2 p' 



dx 2 dy 2 



) 



exactly the same equation that we obtained for the vibrating membrane. 

Remark: The notion of linearization is extremely powerful because it enables 
us to derive information on the behavior of solutions to the nonlinear Eulcr 
equations, which are extremely difficult to solve except for under very special 
circumstances. 

The Euler equations for a perfect gas and the closely related Navier-Stokes 
equations for an incompressible fluid such as water form basic models for fluid 
mechanics. In the case of incompressible fluids, the density is constant, so no 
equation of state is assumed. To allow for viscosity, one adds an additional term 
to the expression (5.17) for the force acting on a fluid clement: 



Here the Laplace operator is applied componentwise and v is a constant, called 
the viscosity of the fluid. The equations used by Navier and Stokes to model an 
incompressible viscous fluid (with p constant) are then 



It is remarkable that these equations, so easily expressed, are so difficult to 
solve. Indeed, the Navier-Stokes equations form the basis for one of the seven 
Millenium Prize Problems, singled out by the Clay Mathematics Institute as 
central problems for mathematics at the turn of the century. If you can show 
that under reasonable initial conditions, the Navier-Stokes equations possess a 
unique well-behaved solution, you may be able to win one million dollars. To 
find more details on the prize offered for a solution, you can consult the web 
address: http: / / www.claymath.org/ millennium/ 

Exercise: 

5.4.1. Show that if the tension and density of a membrane are given by vari- 
able functions T(x,y) and p(x,y) respectively, then the motion of the string is 
governed by the equation 



Force = ^(Av)(x, y, z, 



t)-Vp(x(t),y{t),z{t),t). 



V-v = 0, 



dv 

~dt 



+ (v • V)v = vAw Vp. 

P 



d 2 u i 



V -{T{x,y)Vu). 



dt 2 p(x, y) 
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5.5 Initial value problems for wave equations 

The most natural initial value problem for the wave equation is the following: 

Let I? be a bounded region in the (x, y)-plane which is bounded by a piece- 
wise smooth curve dD, and let hi, hi : D —> R be continuous functions. Find a 
function u(x, y, t) such that 

1. u satisfies the wave equation 

d 2 u 2 fd 2 u d 2 u\ , ro _ 

c 2 • (5-21) 



dt 2 \ dx 2 dy 2 

2. u satisfies the "Dirichlet boundary condition" u(x,y,t) — 0, for (x,y) e 
8D. 

3. u satisfies the initial condition 

u{x,y,0) = hi(x,y), —(x,y,0) = h 2 {x,y). 

Solution of this initial value problem via separation of variables is very similar 
to the solution of the initial value problem for the heat equation which was 
presented in Section 5.3. 

As before, let us suppose that D — {(x, y) € J? 2 : < x < a, < y < b}, so 
that the Dirichlet boundary condition becomes 

u(0, y, t) — u(a, y, t) = u(x, 0, t) = u(x, b, t) = 0. 

We write 

u(x,y,t) = f(x,y)g(t), 
and substitute into (5.21) to obtain 

32 f a2 



Then we divide by c 2 f(x,y)g(t), 



;9"(t) = — - 



c2 g(t) f(x,y) \dx 2 dy 

Once again, we conclude that both sides must be constant. If we let A denote 
the constant, we obtain 



c 2 g{t) f(x,y) \dx 2 dy 2 

This separates into an ordinary differential equation 

g"(t) = c 2 Xg(t), (5.22) 
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and the Helmholtz equation 

The Dirichlct boundary condition becomes 

f(0,y) = f(a,y) = f(x,0) = f(x,b)=0. 

The Helmholtz equation is solved in exactly the same way as before, the 
only nontrivial solutions being 

f mn (x, y) = sm(mirx/a) s'm(niry/b), with X mn = -(mir/a) 2 - (nir/b) 2 . 

The corresponding solution to (5.22) is 

g(t) = Acos(u mn t) + Bsm(u mn t), where u> mn = c^/-X mn . 

Thus we obtain a product solution to the wave equation with Dirichlet boundary 
conditions: 

u(x,y,t) — sm(mirx/a) sin(rt7ry /b)[Acos(ui mn t) + B sm(u> mn t)]. 

The general solution to to the wave equation with Dirichlet boundary conditions 
is a superposition of these product solutions, 

oo 

u(x, y,t) = ^2 sm(rmrx/ a) s'm(niry/b) [A mn cos{u mn t) + B mn sm(u mn t)] . 

m/n—l 

The constants A mn and B mn are determined from the initial conditions. 

The initial value problem considered in this section could represent the mo- 
tion of a vibrating membrane. Just like in the case of the vibrating string, the 
motion of the membrane is a superposition of infinitely many modes, the mode 
corresponding to the pair (m, n) oscillating with frequency ui mn /2n. The lowest 
frequency of vibration or fundamental frequency is 



27r~27rVw \b) ~ 2 



1\ 2 /r 2 

+ 



Exercises: 

5.5.1. Solve the following initial value problem for a vibrating square membrane: 
Find u(x, y,t) 7 0<x<ir,0<y<ir, such that 

d 2 u d 2 u d 2 u 
dt 2 dx 2 dy 2 ' 

u(x, 0, t) = u(x, 7r, t) = u(0, y, t) = u(ir, y, t) = 0, 
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Ou 

u(x,y,0) = 3sinxsiny + 7sin2xsiny, — (x,y, 0) = 0. 

5.5.2. Solve the following initial value problem for a vibrating square membrane: 
Find u(x, y,t) 7 0<x<ir,0<y<ir, such that 

d 2 u (d 2 u d 2 u \ 

W~ \dx^ + d^)' 

u{x, 0, t) = u(x, 7r, t) = u(0, y, t) = u(ir, y, t) = 0, 
du 

u(x,y,0) = 0, — (x,y, 0) = 2 sin x sin y+ 13 sin 2a; sin y. 

5.5.3. Solve the following initial value problem for a vibrating square membrane: 
Find u(x, y,t) 7 0<x<ir,0<y<ir, such that 

8 2 u_d 2 u <Pu 
dt 2 dx 2 dy 2 ' 

u(x, 0, t) — u(x, 7r, t) = u(0, y, t) = u(ir, y, t) = 0, 

du , 



w(a;,y,0) =p(x)q(y), —(x,y,0) = 0, 



where 



\x, for < x < it 12, , . ( y, for < y < tt/2, 

I 7r — x, tor 7T/2 < x < 7r, I 7T — y, lor 7T/2 < y < TT. 

5.6 The Laplace operator in polar coordinates 

In order to solve the heat equation over a circular plate, or to solve the wave 
equation for a vibrating circular drum, we need to express the Laplace opera- 
tor in polar coordinates (r, 0) . These coordinates are related to the standard 
Euclidean coordinates by the formulae 

x = r cos 9, y = r sin 9. 

The tool we need to express the Laplace operator 

dx 2 dy 2 

in terms of polar coordinates is just the chain rule, which we studied earlier in 
the course. Although straigtforward, the calculation is somewhat lengthy. Since 

dx dy d 

— — —(rcos9) = cos 9, — — —(rs'm9) = s'm9, 

or or or or 
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it follows immediately from the chain rule that 
du du dx 



^ dudy _ ^ du ^ du 
dr dx dr dy dr dx dy 



Similarly, since 



d 



- = -<r a »0) = -rsm0 t - = ^(rsinfl) 



rcosf 



it follows that 



du du dx du dy , . . du , . 

de = d-xde + dyde = { - rsm9) d^ + {rcos9) 



du 
dy' 



We can write the results as operator equations, 

Wr = {C ° Sd) lL + (8in0) ^' I = + ^ rC ° s9 ^y 

For the second derivatives, we find that 



d 2 u d du 
dr 2 dr dr 



d_ 
dr 



n du . n du 
cos ft- — h sint/— 
dx dy 



= cos 9 



d f du 



dr \ dx 



+ sin 9 



d f du 



dr \dy 



Similarly, 



2 „a 2 u „ n . . d 2 u . 2 n d 2 u 
= cos 0— — + 2 cos 6* sin 6*— — - + sin V— — . 
dx dxdy dy z 



d 2 u d du d 
W = WW = d9 



. n du n du 
—rsmo- — h rcosu — 
dx dy 



d f du 



-rsini 



d9 V dx 



d f du 



r cos ( 



,aw 



,du 



d9 \dy 

2 . 9 r,d 2 U „ 2 r, ■ n 9 2 U 

r sm 9— — 2r cos 9 sm 9 — — - 
ax z array 



rcosf — rsmtf — — 

ax aw 



which yields 

1 a 2 u 



• 2 „<9 2 w „ . . . d 2 u 
r l dv z dx z dxdy 



cos 2 9 



■cos 2 9 1 — , 
dy 1 

d 2 u 1 du 



du 
dr ' 



dy 2 



r dr 



Adding these results together, we obtain 



d 2 u 1 d 2 u d 2 u d 2 u 1 du 
dr 2 r 2 d9 2 dx 2 dy 2 r dr ' 



or equivalently, 



Au _ d 2 u 1 d 2 u 1 du 
dr 2 r 2 d9 2 r dr 
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Finally, we can write this result in the form 

1 d ( du\ 1 8 2 u 

This formula for the Laplace operator, together with the theory of Fourier 
series, allows us to solve the Dirichlet problem for Laplace's equation in a disk. 
Indeed, we can now formulate the Dirichlet problem as follows: Find u(r, 0), for 
< r < 1 and 6 mathbbR, such that 

1. u satisfies Laplace's equation, 

r dr \ dr ) r 2 dO 2 ' 

2. u satisfies the periodicity condition u(r, 9 + 2n) = u(r, 9), 

3. u is well-behaved near r — 0, 

4. u satisfies the boundary condition u(l,6) = h{9), where h(9) is a given 
well-behaved function satisfying the periodicity condition h{9-\-2ir) = h{9). 

The first three of these conditions are homogeneous linear. To treat these 
conditions via the method of separation of variables, we set 

u(r, 9) = R(r)0(9), where 6(0 + 2tt) = 6(0). 

Substitution into (5.25) yields 

1 d ( dR\ Rd 2 _ () 
r dr \ dr J r 2 d9 2 



We multiply through by r 



2 



d ( dR\^ d 2 B 
r— r— 6 + R- 



dr V dr J d9 2 
and divide by — RQ to obtain 

r d ( dR\ 1 d 2 Q 



Rdr \ dr J 6 d9 2 ' 

The left-hand side of this equation does not depend on while the right-hand 
side does not depend on r. Thus neither side can depend on cither or r, and 
hence both sides must be constant: 

r d f r dR\ _ 1_^6 _ 



Rdr \ dr J 6 d9 2 

Thus the partial differential equation divides into two ordinary differential 
equations 

^ = A6, 6(0 + 270 = 6(0), 
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We have seen the first of these equations before when we studied heat flow in a 
circular wire, and we recognize that with the periodic boundary conditions, the 
only nontrivial solutions are 

a = o, e = f, 

where ao is a constant, and 

A = — n 2 , O = a n cos n9 + b n sin n9, 

where a n and b n are constants, for n a positive integer. Substitution into the 
second equation yields 

d ( dR\ 2 
r — r — - n 2 R = 0. 
dr \ dr J 

If n — 0, the equation for R becomes 

d ( r dR\ _ 
dr \ dr J 

which is easily solved to yield 

R(r) = A + B log r, 

where A and B are constants of integration. In order for this equation to be well- 
behaved as r — > we must have B = 0, and the solution uo(r,6) to Laplace's 
equation in this case is constant. 

When n^O, the equation for R is a Cauchy-Eulcr cquidimensional equation 
and we can find a nontrivial solution by setting 

R(r) = r m . 

Then the equation becomes 

r4~ ( r4- (r m )) - n 2 r m 7 
dr \ dr J 

and carrying out the differentiation on the left-hand side yields the characteristic 
equation 

m 2 — n 2 = 0, 

which has the solutions m = ±n. In this case, the solution is 

R(r) = Ar n + Br- n . 

Once again, in order for this equation to be well-behaved as r — > we must 
have B = 0, so R(r) is a constant multiple of r n , and 



u n (r,8) — a n r n cos n9 + b n r n sin n6. 
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The general solution to (5.25) which is well-behaved at r = and satisfies 
the periodicity condition u(r, 9 + 2tt) = u(r, 9) is therefore 

oo 

u(r, 9) = — + ^ (a n r n cos n6 + b n r n sin nO) , 

71=1 

where a , a 1; . . . , b\, . . . are constants. To determine these constants we must 
apply the boundary condition: 

oo 

h(9) = u(l,9) = + ^(a n cosnd + b n sinn9). 
2 „=i 

We conclude that the constants a , Oi, . . . , b\, . . . are simply the Fourier coeffi- 
cients of h. 

Exercises: 

5.6.1. Solve the following boundary value problem for Laplace's equation in a 
disk: Find u(r, 9),0 < r < 1, such that 

1 d ( du\ 1 d 2 u 

r dr \ dr ) r 2 d9 2 ' 

u(r, 9 + 2ir) = u(r, 9), u well-behaved near r = 

and 

u(l,6) = h(6), where h{9) = 1 + cos6» - 2sin6» + 4 cos 2(9. 

5.6.2. Solve the following boundary value problem for Laplace's equation in a 
disk: Find u(r,0),0 < r < 1, such that 



ld_ ( du 
r dr \ dr / 

u(r, 9 + 2n) = u(r, 9), u well-behaved near r = 



r — I _|_ _L = o, 
dr V dr J r 2 d9 2 



and 

u(l,6) = h(6), 
where h{9) is the periodic function such that 

h{6) = \0\, for -7T < 6> < tt. 

5.6.3. Solve the following boundary value problem for Laplace's equation in an 
annular region: Find u(r,9), 1 < r < 2, such that 

i a / au\ ia 2 «_ o 

r dr \ dr J r 2 d9 2 
u(r,9 + 2ir) =u(r,0), u(l,6) = 1 + 3 cos 9 - sin 9 + cos 29 

and 

u(2,6) = 2cos6» + 4cos26». 
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5.7 Eigenvalues of the Laplace operator 

We would now like to consider the heat equation for a room whose shape is given 
by a well-behaved but otherwise arbitrary bounded region D in the {x, i/)-plane, 
the boundary dD being a well-behaved curve. We would also like to consider 
the wave equation for a vibrating drum in the shape of such a region D. Both 
cases quickly lead to the eigenvalue problem for the Laplace operator 

A- _^- + ^_ 
dx 2 dy 2 

Let's start with the heat equation 

^ = c 2 A U , (5.26) 
with the Dirichlct boundary condition, 

u(x,y,t)=0 for (x 1 y)€dD. 
To solve this equation, we apply separation of variables as before, setting 

u(x,y,t) = f{x,y)g(t), 
and substitute into (5.26) to obtain 

f(x,y)g'(t) = c 2 (Af)(x,y)g(t). 
Then we divide by c 2 f(x, y)g(t), 

' g'(t) = -^—<Af)(x,y). 



c 2 g{t) f(x,y) 

The left-hand side of this equation does not depend on x or y while the right- 
hand side does not depend on t. Hence both sides equal a constant A, and we 
obtain 

^' (t) = A = 7(^) (a/)feri - 

This separates into 

g'(t) = c 2 \g{t) and (A/)(x, y) = Xf(x, y) 
in which / is subject to the boundary condition, 

f(x,y) = for (x,y)edD. 
The same method can be used to treat the wave equation 

^ = c 2 A., (5.27) 



141 



with the Dirichlet boundary condition, 

u(x,y, t)=0 for (x,y) £ dD. 
This time, substitution of u(x,y,t) — f(x,y)g(t) into (5.27) yields 

f(x,y)g"(t)=c 2 (Af)(x,y)g(t), or ^L- g '' {t ) = j^-^Af^y). 

Once again, both sides must equal a constant A, and we obtain 

This separates into 

g"(t) = Xg(t) and (Af)(x, y) = Xf(x, y), 

in which / is once again subject to the boundary condition, 

f(x,y) = for (x,y)edD. 

In both cases, we must solve the "eigenvalue problem" for the Laplace op- 
erator with Dirichlet boundary conditions. If A is a real number, let 

W x = {smooth functions f(x, y) : Af = A/, f\dD = 0}. 

We say that A is an eigenvalue of the Laplace operator A on D if W\ ^ 0. 
Nonzero elements of W\ are called eigenfunctions and W\ itself is called the 
eigenspace for eigenvalue A. The dimension of W\ is called the multiplicity of 
the eigenvalue A. The eigenvalue problem consist of finding the eigenvalues A, 
and a basis for each nonzero eigenspace. 

Once we have solved the eigenvalue problem for a given region D in the 
(x, y)-planc, it is easy to solve the initial value problem for the heat equation 
or the wave equation on this region. To do so requires only that we substitute 
the values of A into the equations for g. In the case of the heat equation, 

g'(t) = c 2 \g(t) g(t) = (constant)e c2At , 

while in the case of the wave equation, 

g"(t) = c 2 Xg(t) g(t) = (constant) sin(c\/ z A(t - t )). 

In the second case, the eigenvalues determine the frequencies of a vibrating 
drum which has the shape of D. 

Theorem. All of the eigenvalues of A are negative, and each eigenvalue has 
finite multiplicity. The eigenvalues can be arranged in a sequence 

> Ai > A 2 > • • • > A„ > • • • , 
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with X n — » — oo. Every well-behaved function can be represented as a convergent 
sum of eigcn functions. 2 

Although the theorem is reassuring, it is usually quite difficult to determine the 
eigenvalues Ai, A 2 , • • • , A„, • • • explicitly for a given region in the plane. 

There are two important cases in which the eigenvalues can be explicitly 
determined. The first is the case where D is the rectangle, 

D = {(x, y) : < x < a, < y < b}. 

We saw how to solve the eigenvalue problem 

+ = A /' /(*.») = « fOT (x,y)€dD, 

when we discussed the heat and wave equations for a rectangular region. The 
nontrivial solutions are 

A - f™ 1 } 2 (^Y f (r - sin ( l:mc \ sin ( im V\ 

where m and n are positive integers. 

A second case in which the eigenvalue problem can be solved explicitly is 
that of the disk, 

D = {(x,y) :x 2 + y 2 < a 2 }, 
where a is a positive number, as we will see in the next section. 3 
Exercises: 

5.7.1. Show that among all rectangular vibrating membranes of area one, the 
square has the lowest fundamental frequency of vibration by minimizing the 
function 

'<-')-□'+(!)' 

subject to the constraints ab = 1, a, b > 0. Hint: One can eliminate b by setting 
a = 1/b and then find the minimum of the function 

9(a) = (^) +Oa) 2 , 

when a > 0. 

5.7.2. Let D be a finite region in the (a;,y)-plane bounded by a smooth curve 
3D. Suppose that the eigenvalues for the Laplace operator A with Dirichlct 
boundary conditions on D arc Ai, A 2 , . . . , where 

> Ai > A 2 > • • • , 

2 Note the similarity between the statement of this theorem and the statement of the 
theorem presented in Section 4.7. In fact, the techniques used to prove the two theorems are 
also quite similar. 

3 A few more cases are presented in advanced texts, such as Courant and Hilbcrt, Methods 
of mathematical physics I , New York, Intcrscicnce, 1953. See Chapter V, §16.3. 
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each eigenvalue having multiplicity one. Suppose that <p n (x,y) is a nonzero 
cigenfunction for eigenvalue A„. Show that the general solution to the heat 
equation 

du 

di= Au 

with Dirichlct boundary conditions (u vanishes on dD) is 

oo 

u{x,t) = ^2b n <j) n (x,y)e > 



where the b n 's are arbitrary constants. 

5.7.3. Let D be a finite region in the (x, y)-plane as in the preceding problem. 
Suppose that the eigenvalues for the Laplace operator A on I? are Ai, A 2 , . . . , 
once again. Show that the general solution to the wave equation 

d 2 u A 

together with Dirichlet boundary conditions (u vanishes on dD) and the initial 
condition 

|W,0) = 0, 

is 

oo 

u ( x , t) = h n<t>n{x, V) COs( V / -A„t), 

n=l 

where the 6„'s are arbitrary constants. 

5.8 Eigenvalues of the disk 

To calculate the eigenvalues of the disk, it is convenient to utilize polar coordi- 
nates r, 9 in terms of which the Laplace operator is 



1 d ( df\ 1 d 2 f 



r dr \ dr I r 2 d6 2 



= A/, f\8D = 0. (5.28) 



Once again, we use separation of variables and look for product solutions of the 
form 

f(r,6) = R(r)Q(6), where R(a) = 0, 0(0 + 2tt) = 0(0). 
Substitution into (5.28) yields 

r dr \ dr ) r 2 dQ 2 
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We multiply through by r 



2 



d ( dR\^ d 2 



and divide by RQ to obtain 



r d ( dR\ ~ 1 <PG 
r— — — \r 



Rdr V dr J 6 d6 2 ' 

Note that in this last equation, the left-hand side does not depend on 9 while 
the right-hand side does not depend on r. Thus neither side can depend on 
either 9 or r, and hence both sides must be constant: 

r d ( dR\ - 1 d 2 @ 

~RTr\ r ^) +Xr =e^ = /i - 

Thus in the manner now familiar, the partial differential equation divides 
into two ordinary differential equations 

^ = m6, e(9 + 2n) = e(9), 

r d / dR\ _ Xr2R= _ R( a )=0. 
dr \ dr J 

Once again, the only nontrivial solutions are 

m = o, e = f 

and 

fj,— —n 2 , Q — a n cosn9 + b n sinn9, 
for n a positive integer. Substitution into the second equation yields 

d fdR\ +{ _ Xr 2_ n 2 )R = 

dr \ dr J 

Let x = \/—\r so that x 2 = — Xr 2 . Then since 



dr c?x' 



our differential equation becomes 

d/dR\ +3 2_ n2 
dx \ dx J 

where R vanishes when x = ^f—\a. If we replace R by y, this becomes 

+ (x 2 - n 2 )y = 0, y(V^Xa) = 0. (5.30) 



d ( dy 

x 



dx V dx 
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The differential equation appearing in (5.30) is Bessel's equation, which we 
studied in Section 1.4. 

Recall that in Section 1.4, we found that for each choice of n, Bessel's equa- 
tion has a one-dimensional space of well-behaved solutions, which are constant 
multiples of the Bessel function of the first kind J n (x). Here is an important 
fact regarding these Bessel functions: 

Theorem. For each nonnegative integer n, J n (x) has infinitely many positive 
zeros. 

Graphs of the functions Jo(x) and J\(x) suggest that this theorem might well be 
true, but it takes some effort to prove rigorously. For completeness, we sketch 
the proof for the case of Jq(x) at the end of the section. 

The zeros of the Bessel functions are used to determine the eigenvalues of 
the Laplace operator on the disk. To see how, note first that the boundary 
condition, 

y(>/=Aa) = 0, 

requires that \/— Xa be one of the zeros of J n (x). Let a n:k denote the fc-th 
positive root of the equation J n (x) = 0. Then 

a 2 

V-Xa = a n ^ k => A= 7 y r , 

and 



a 2 



R{r) = J n (a n ,kr/a) 

will be a solution to (5.29) vanishing at r = a. Hence, in the case where n = 0, 



and 

/o,fc(r,0) = J (a , k r/a) 
will be a solution to the eigenvalue problem 



r dr V drj^r 2 dd 2 ^ ^ 



Similarly, in the case of general n, 

\ - Q ".fc 

■*-n,k 9 ; 

and 

fn,k{r, 9) = J n (a nik r/a)cos(n8) or g n>k = J n (a rltk r/a) sm(n0), 
for n = 1,2, ... , will be solutions to this eigenvalue problem. 
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Figure 5.2: Graph of the function fo,i(r, 0). 



Each of these eigenfunctions corresponds to a mode of oscillation for the 
circular vibrating membrane. Suppose, for example, that a = 1. Then 

a QA = 2.40483 => A ,i = -5.7832, 

a ,2 = 5.52008 =*> A , 2 = -30.471, 

a ,3 = 8.65373 A , 3 = -74.887, 

a ,4 = 11.7195 => A ,3 = -139.039, 

and so forth. The mode of oscillation corresponding to the function /o,i will 
vibrate with frequency 

\Z~Vi _ Qui 
2tt ~ 2?r 



= .38274, 



while the mode of oscillation corresponding to the function /o,2 will vibrate with 
frequency 

\J — A ,2 _ «0,2 



2?r 2tt 



= .87855. 



Similarly, 



a x> \ = 3.83171 Ai,i = -14.682, 

a.2,1 = 5.13562 => A 2 ,i = -26.3746, 

and hence the mode of oscillation corresponding to the functions and g\ \ 
will vibrate with frequency 



2tt 2tt 



= .60984, 



147 



1 

Figure 5.3: Graph of the function fo,2(r,0)- 



while the mode of oscillation corresponding to the functions /2.1 and 172,1 will 
vibrate with frequency 

^Ihl = = .81736. 
2tt 2tt 

A general vibration of the vibrating membrane will be a superposition of 
these modes of oscillation. The fact that the frequencies of oscillation of a 
circular drum are not integral multiples of a single fundamental frequency (as 
in the case of the violin string) limits the extent to which a circular drum can 
be tuned to a specific tone. 

Proof that Jq(x) has infinitely many positive zeros: First, we make a change 
of variables x — e z and note that as z ranges over the real numbers, the corre- 
sponding variable x ranges over all the positive real numbers. Since 

did d 1 d d 
dx = e dz, — — — — and hence x— — = e — = — . 

dx e z dz dx e z dz dz 



Thus Bessel's equation (5.30) in the case where n = becomes 

+ e 2z y = 0. (5.31) 



dz 2 



Suppose that z > 1 and y(z) is a solution to (5.31) with y(z ) ^ 0. We 
claim that y(z) must change sign at some point between zq and zq + n. 
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Figure 5.4: Graph of the function /i,i(r, ff). 



Assume first that y(zo) > 0. Let f(z) = sm(z — zq). Since f"(z) = —f(z) 
and y(z) satisfies (5.31), 

± [y(z)f'(z) - y'(z)f(z)} = y(z)f"(z) - y"{z)f{z) 

= -y(z)f(z) + e 2z y(z)f(z) = (e 2 * - l)y(z)f(z). 

Note that f(z) > for z between zq and zq + it. If also y(z) > for z between 
zq and zq + 7r, then y{z)f{z) > and 

< / (e 2z - l)y(z)f(z)dz = [y{z)f{z) - y> \z) f (z)] z ° + * 

= v{zq + 7r)/'(zo + 7r) - y(z )f(z ) = -y(z + 7r) - y(z ) < 0, 

a contradiction. Hence our assumption that y(z) be postive for z between zq 
and zq + 7r must be incorrect, y(z) must change sign at some z between Zq and 
zq + 7T, and hence y{z) must be zero at some point in this interval. 

If y{zo) < 0, just apply the same argument to —y. In either case, we conclude 
that y(z) must be zero at some point in any interval to the right of z = 1 of 
length at least 7r. It follows that the solution to (5.31) must have infinitely many 
zeros in the region z > 1, and Jo (a:) must have infinitely many positive zeros, 
as claimed. 

The fact that J n (x) has infinitely many positive zeros could be proven in a 
similar fashion. 
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Exercises: 

5.8.1. Solve the following initial value problem for the heat equation in a disk: 
Find u(r, B,t),0 < r < 1, such that 

du _ 1 d ( du\ 1 d 2 u 
dt r dr \ dr J r 2 d9 2 ' 

u(r, 9 + 2ir, t) — u(r, 6,t), u well-behaved near r = 0, 

u{i,e,t) = o, 

u(r,9,0) = Jo{ao,ir) + 3Jo(a Q ^r) - 2Ji{ot\ t \r) sin 6* + AJ 2 {a 2 ,ir) cos 29. 

5.8.2. Solve the following initial value problem for the vibrating circular mem- 
brane: Find u(r, 6, t), < r < 1, such that 

d 2 u _ 1 d f du\ 1 d 2 u 
dt 2 r dr \ dr ) r 2 89 2 ' 

u(r, 9 + 2-k, t) — u(r, 6,t), u well-behaved near r — 0, 

u(l,M)=0, ^(r,0,O) = O, 
u(r, 0,0) = Jo(ao,i r ) + 3Jo(ao,2T") - 2Ji(ai,ir) sin 6* + 4J 2 (o!2,ir) cos 20. 

5.8.3. (For students with access to Mathematica) a. Run the following Mathe- 
matica program to sketch the Bessel function Jq(x): 

n=0; Plot [ Bessel J [n,x] , {x,0,15}] 

b. From the graph it is clear that the first root of the equation Jq(x) — is 
near 2. Run the following Mathematica program to find the first root ao,i of 
the Bessel function Jo (a;): 

n=0; FindRoot[ Bessel J [n,x] == ,{x,2}] 

Find the next two nonzero roots of the Bessel function Jq(x). 

c. Modify the programs to sketch the Bessel functions Ji(x),... ,J$(x), and 
determine the first three nonzero roots of each of these Bessel functions. 

5.8.4. Which has a lower fundamental frequency of vibration, a square drum or 
a circular drum of the same area? 

5.9 Fourier analysis for the circular vibrating 
membrane* 

To finish up the solution to the initial value problem for arbitrary initial displace- 
ments of the vibrating membrane, we need to develop a theory of generalized 
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Fourier series which gives a decomposition into cigenfunctions which solve the 
eigenvalue problem. 

Such a theory can be developed for an arbitrary bounded region D in the 
(x, y)-plane bounded by a smooth closed curve dD. Let V denote the space of 
smooth functions / : D — > R whose restrictions to dD vanish. It is important 
to observe that V is a "vector space" : the sum of two elements of V again lies 
in V and the product of an clement of V with a constant again lies in V. 

We define an inner product (, ) on V by 



Lemma. With respect to this inner product, eigenfunctions corresponding to 
distinct eigenvalues are perpendicular; if f and g are smooth functions vanishing 
on dD such that 



then cither X — fi or (/, g) = 0. 

The proof of this lemma is a nice application of Green's theorem. Indeed, it 
follows from Green's theorem that 




A/ = A/, Ag = fig, 



(5.32) 




Hence if the restriction of f to dD vanishes 




Similarly if the restriction of g to dD vanishes, 




Thus if / and g lie in V, 



(/, A 5 ) = - / / V/ • Vgdxdy - (g, A/) 



J Jd 



In particular, if (5.32) holds, then 



n(f,g) = (f,Ag) = (Af,g) = \(f,g) 



and hence 



(A-/i)</,0>=O. 



151 



It follows immediately that either A — //, = or (/, g) = 0, and the lemma is 
proven. 

Now let us focus on the special case in which D is a circular disk. Recall the 
problem that we want to solve, in terms of polar coordinates: Find 

u(r,6,t),0 < r < 1, 

so that 

d 2 u 1 d ( du\ 1 d 2 u 



dt 2 r dr \^ dr ) r 2 d9 2 ' 
u(r, 9 + 2ir, t) — u(r, 9,t), u well-behaved near r = 0, 

u(i,e,t) = o, 

u(r,6,0) = h(r,6), — (r,fl,0)=0. 

It follows from the argument in the preceding section that the general solu- 
tion to the homogeneous linear part of the problem can be expressed in terms 
of the cigenfunctions 

fn,k{r, 9) = J n {a n ,kr/a) cos(n#), g nik = J n (a n , k r/a) sm(n6). 

Indeed, the general solution must be a superposition of the products of these 
eigenfunctions with periodic functions of t of frequencies y/— X ny k/2n. Thus this 
general solution is of the form 



u ( r , 9 ^) = Y1 a o,fc5o,fcO, 9) cos(a ,fet) 



fc=i 



+ ^^K,fc/n,fc(r,0) + & n ,fe.9 n ,fe(r, 9)] cos(a n . k t). (5.33) 

n=l fe=l 

We need to determine the a n ,kS and the b n ,kS so that the inhomogencous initial 
condition u(r, 9,0) = h(r,9) is also satisfied. If we set t = in (5.33), this 
condition becomes 

oo oo oo 

^2 ao,kfo,k + ^2 ^2[a n ,kfn,k + b n , k g n ,k} = h. (5.34) 

fc=l n=lfe=l 

Thus we need to express an arbitrary initial displacement h(r, 9) as a superpo- 
sition of these eigenfunctions. 

It is here that the inner product (, ) on V comes to the rescue. The lemma 
implies that cigenfunctions corresponding to distinct eigenvalues are perpendic- 
ular. This implies, for example, that 

(fo,jJo,k) = 0, unless j = k. 
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Similarly, for arbitrary n > 1, 

{fn,j,fn,k) = (fn,j,g n ,k) = {9n,j,9n,k) = 0, unless j = k. 

The lemma docs not immediately imply that 

(fn,j ? 9n,j) 0. 

To obtain this relation, recall that in terms of polar coordinates, the inner 
product (,) on V takes the form 



(f,9) 



1 r /-2ir 



f(r,e)g(r,6)rd6 



dr. 



If we perform integration with respect to 6 first, we can conclude from the 
familiar integral formula 



2tt 

sin n9 cos n9d9 = 



/>Z7T 
/ ^ 

Jo 



that 

(/nji 9n,j) = 0, 

as desired. 

Moreover, using the integral formulae 



I 



2tt 

cos n9 cos mOdO 



2tt 

sin n9 sin mOdO 



| 7r, for m = n, 

I 0, for m ^= n, 

7r, for m = n, 

0, for m^= n, 



L 



2tt 

sin n6 cos mOdO = 0, 



we can check that 

(fnj:fm,k) — {9n,j:9m,k) — (fn,ji9m,k) 0: 

unless m = n. It then follows from the lemma that 

(/„,,-, fm,k) - (g n j,9m,k) = (fn,j,9m,k) = 0, 

unless j = k and m — n. 

From these relations, it is not difficult to construct formulae for the coeffi- 
cients of the generalized Fourier series of a given function h(r, 9). Indeed, if we 
take the inner product of equation (5.34) with u n ^k, we obtain 

Q>n,k {fn,k ? fn.k) (^-?/n,fc)i 
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or equivalcntly 



(h,f n ,k) 



(5.35) 




Similarly, 



b n ,k 



(h,g n ,k) 



(5.36) 



(<7n,fcj <7n,fc) 



Thus wc finally see that the solution to our initial value problem is 



oc 



u(r, 0,t) = y~] a ,fe/o,fc(r, 6) cos(a ,fc£) 



oo 



oo 



n=l fe=l 



(r, 0) + b n , k g n ,k{r, 9)] cos(a n , fe i), 



where the coefficients are determined by the integral formulae (5.35) and (5.36). 
Exercises: 

5.9.1. (For students with access to Mathematica) Suppose that a circular vi- 
brating drum is given the initial displacement from equilibrium described by 
the function 



for < r < 1 . In order to find the solution to the wave equation with this initial 
condition, we need to expand h(r, 9) in terms of eigenfunctions of A. Because 
of axial symmetry, we can write 



a. Use the following Mathematica program to determine the coefficient a 01 in 
this expansion: 

d = NIntegrate[r (BesselJ [0,2. 40483 r] ) A2 , {r , , 1}] ; 

n = NIntegrate[r (BesselJ [0,2. 40483 r] ) .1 (1-r) ,{r,0,l}] ; 

aOl = n/d 

b. Use a modification of the program to determine the coefficients ao,2, ao,3, 
and a ,4. 

c. Determine the first four terms in the solution to the initial value problem: 



h(r,9) = .1(1 -r) 



h(r, 9) = a a i J (a ^r) + a ^Jo(ao,2r) + . . . . 



Find 



u(r,9,t),0 < r < 1, 



so that 
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u(r, 9 + 2ir, t) — u(r, 6,t), u well-behaved near r = 0, 

u(i,e,t) = o, 

OIL 

u(r,fl,0) = .l(l-r), — M,0) = 0. 

5.9.2. (For students with access to Mathematica) The following program will 
sketch graphs of the successive positions of a membrane, vibrating with a su- 
perposition of several of the lowest modes of vibration: 

« Graphics 'ParametricPlot3D' 

aOl = .2; a02 = .1; a03 = .2; a04 = .1; all = .3; bll = .2; 

vibmem = Table [ 

CylindricalPlot3D [ 

aOl BesselJ[0,2.4 r] Cos[2.4 t] 
+ a02 BesselJ[0,5.52 r] Cos [5. 52 t] 
+ a03 BesselJ[0,8.65 r] Cos [8. 65 t] 
+ a04 BesselJ[0,11.8 r] Cos [11. 8 t] 
+ all BesselJ[l,3.83 r] Cos [u] Cos[3.83 t] 
+ bll BesselJ[l,3.83 r] Sin[u] Cos[3.83 t] , 
{r,0,l}, {u,0,2 Pi}, PlotRange -> {-.5,. 5} 
], {t,0,l.,.l} 

] 

a. Execute this program and then animate the sequence of graphics by running 
"Animate selected graphics," from the Graphics menu. (Warning! Generating 
the sequence of graphics takes lots of memory. If you have insufficient memory 
for the animation, try inserting PlotPoints -> 8 within the CylindricalPlot 
statement.) 

b. Replace the values of aOl, . . . with the Fourier coefficients obtained in the 
preceding problem, execute the resulting program to generate a sequence of 
graphics, and animate the sequence as before. 
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Appendix A 



Using Mathematica to solve 
differential equations 

In solving differential equations, it is sometimes necessary to do calculations 
which would be prohibitively difficult to do by hand. Fortunately computers 
can do the calculations for us, if they are equiped with suitable software, such 
as Maple, Matlab, or Mathematica. This appendix is intended to give a brief 
introduction to the use of Mathematica for doing such calculations. 

Most computers which are set up to use Mathematica contain an on-line 
tutorial, "Tour of Mathematica," which can be used to give a brief hands-on 
introduction to Mathematica. Using this tour, it may be possible for many 
students to learn Mathematica without referring to lengthy technical manuals. 
However, there is a standard reference on the use of Mathematica, which can 
be used to answer questions if necessary. It is The Mathematica Book, by 
Stephen Wolfram, Fourth edition, Wolfram Media and Cambridge University 
Press, 1999. 

We give here a very brief introduction to the use of Mathematica. After 
launching Mathematica, you can use it as a "more powerful than usual graphing 
calculator." For example, if you type in 

(11 - 5)/3 

the computer will perform the subtraction and division, and respond with 
Dut[l] = 2 

The notation for multiplication is *, so if you type in 

2 * (7 + 4) 

the computer will respond with 
Out [2] = 22 
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You can also use a space instead of * to denote multiplication, so if you input 

2 8 

the computer will respond with 
Out [3] = 16 

The computer can do exact arithmetic with integers or rational numbers. 
For example, since A is the symbol used for raising to a power, if you type in 

2A150 

the computer will calculate the 150th power of two and give the exact result: 

Out [4] = 1427247692705959881058285969449495136382746624 

On the other hand, the symbol N tells the computer to use an approximation 
with a fixed number of digits of precision, so entering 

N[2A150] 

will give an approximation to the exact answer, expressed in scientific notation: 

Out [5] = 1.42725 10 45 

Real numbers which are not rational, such as n, have infinite decimal expan- 
sions which cannot be stored within a computer. However, we can approximate 
a real number by a finite decimal expansion with a given number of digits of pre- 
cision. For example, since Mathematica uses the name E to denote the number 
e, typing in 

N[E] 

will tell Mathematica to give a rational approximation to e to a standard number 
of digits of precision: 

Out [6] = 2.71828 

In principle, the computer can calculate rational approximations to an arbitrary 
number of digits of precision. Since the number n is represented by Mathematica 
as Pi, typing in 

N[Pi,40] 

will tell Mathematica to give a rational approximation to n to 40 digits of 
precision: 

Out [7] = 3.1415926535897932384626433832795028841972 

The computer can represent not only numbers but functions. Within Math- 
ematica, built-in functions are described by names which begin with capital 
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Figure A.l: Graph of the logarithm function. 

letters. For example, Log denotes the natural base e logarithm function. Thus 
entering 

N[Log[2] ,40] 

will give the logarithm of 2 with 40 digits of precision: 

Out [8] = 0.693147180559945309417232121458176568076 

One can also plot functions with Mathcmatica. For example, to plot the 
logarithm function from the values 1 to 3, one simply inputs 

Plot [Log[t] ,{t,l,3}] 

and Mathematica will automatically produce the plot shown in Figure A.l. 

We can also define functions ourselves, being careful not to capitalize them, 
because capitals are reserved by Mathematica for built-in functions. Thus we 
can define the function 

y(t) = ce kt 

by typing 

y [t_] := c * EA(k * t) ; 

Mathematica will remember the function we have defined until we quit Mathe- 
matica. We must remember the exact syntax of this example (use of the under- 
line character and the colon followed by the equal sign) when defining functions. 
In this example c and k are parameters which can be defined later. Just as in 
the case of functions, we must be careful to represent parameters by lower case 
letters, so they will not be confused with built-in constants. Further entry of 

c = 1; k = 1; N[y[l]] 
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yields the response 
Out [11] = 2.71828 

while entering 

Plot [y[t] ,{t,0,2}] 

will give a plot of the function we have defined, for t in the interval [0, 2]. 
We can use Mathematica to solve matrix differential equations of the form 

f = .4*. (A.1) 

where A is a square matrix with constant entries. 

The first step consists of using Mathematica to find the eigenvalues and 
eigenvectors of A. To see how this works, we must first become familiar with 
the way in which Mathematica represents matrices. Since Mathematica reserves 
upper case letters for descriptions of built-in functions, it is prudent to denote 
the matrix A by lower case a when writing statements in Mathematica. The 
matrix 

A ={-? -\ 

can be entered into Mathematica as a collection of row vectors, 

a = {{-2, 5}, {1,-3}} 
with the computer responding by 
Dut[l] = {{-2,5}, {1,-3}} 

Thus a matrix is thought of as a "vector of vectors." Entering 

MatrixForm [a] 

will cause the computer to give the matrix in the familiar form 

Out [2] =-2 5 

1 -3 

To find the eigenvalues of the matrix A, we simply type 

Eigenvalues [a] 
and the computer will give us the exact eigenvalues 

—5 — -n/21 -5 + V2T 



which have been obtained by using the quadratic formula. Quite often numerical 
approximations arc sufficient, and these can be obtained by typing 

eval = Eigenvalues [N [a] ] 
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the response this time being 

Out [4] = {-4.79129, -0.208712} 

Defining eval to be the eigenvalues of A in this fashion, allows us to refer to 
the eigenvalues of A later by means of the expression eval. 

We can also find the corresponding eigenvectors for the matrix A by typing 

evec = Eigenvectors [N [a] ] 

and the computer responds with 

Out [5] = {{-0.873154, 0.487445}, {0.941409, 0.337267}} 

Putting this together with the eigenvalues gives the general solution to the 
original linear system (A.l) for our choice of the matrix A: 

v(A r f-0-873154A 4 .79129t , „ f0.941409\ 0.20871M 

x(<) - Cl (, 0.487445 ) e +C2 (,0.337267 ) 6 

Mathematica can also be used to find numerical solutions to nonlinear differ- 
ential equations. The following Mathematica programs will use Mathematica's 
differential equation solver (which is called up by the command NDSolve), to 
find a numerical solution to the initial value problem 

dy/dx = y(l-y), y(0) = .1, 

give a table of values for the solution, and graph the solution curve on the 
interval < x < 6. The first step 

sol := NDSolve[{ y> [x] == y[x] (1 - y [x] ) , y[0] == .1 }, 
y, {x,0,6}] 

generates an "interpolation function" which approximates the solution and calls 
it sol, an abbreviation for solution. We can construct a table of values for the 
interpolation function by typing 

Table [Evaluate [y[x] /. sol], {x,0,6,.l}]; 

or graph the interpolation function by typing 

Plot [Evaluate [y [x] /. sol], {x,0,6}] 

This leads to a plot like that shown in Figure A. 2. 

Readers can modify these simple programs to graph solutions to initial value 
problems for quite general differential equations of the canonical form 

*-/(..,). 

All that is needed is to replace the first argument of NDSolve with the differential 
equation one wants to solve, remembering to replace the equal signs with double 
equal signs, as in the example. 
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1 2 3 4 5 6 

Figure A. 2: Solution to y' = y(l — y), y(0) = .1. 




In fact, it is no more difficult to treat initial value problems for higher order 
equations or systems of differential equations. For example, to solve the initial 
value problem 

^ = -xv, ^ = x y-y> *( ) = 2 > y(o)-- 1 - ( A - 2 ) 

one simply types in 

sol := NDSolve[{ x' [t] == - x[t] y[t], y' [t] == x[t] y[t] - y[t], 
x[0] == 2, y[0] == .1 }, 
{x,y}, {t,0,10}] 

Once we have this solution, we can print out a table of values for y by entering 
Table [Evaluate [y[t] /. sol], {t,0,2,.l}] 
We can also plot y(t) as a function of t by entering 
Plot [Evaluate [y [t] /. sol], {t,0,10}] 

Figure A. 3 shows a parametric plot of (x(t),y(t)) which is generated by entering 
ParametricPlot [Evaluate [{ x[t], y[t]} /. sol], {t,0,10}] 
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